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Particle Splitting: A New Method 
for SPH Star Formation 
Simulations 



Abstract 

Particle Splitting is a new algorithm invented for use with self-gravitating SPH codes. It 
is designed to enable numerical simulations to obey the Jeans condition at all times (but 
it could be used in other contexts, to satisfy other conditions which require high resolution 
locally). With particle splitting, all coarse particles in regions of interest, are erased and 
replaced by sets of fine particles, increasing the numerical resolution of the simulations. A 
new algorithm for calculating smoothing lengths was added to our code, to accommodate the 
mixing of different mass particles. Our particle splitting SPH code reproduced the results of 
standard test simulations. 

Simulations of rotating clouds with m=2 density perturbations produce a binary and a bar. 
We confirm that fragmentation of the bar should be attributed to inadequate rcsohition. By 
applying Particle Splitting to such simulations, we reproduce the results of high resolution 
finite difference simulations (Truelove et al. and Klein et al), where bar fragmentation is 
absent. We obtain these results with great computational economy. 

We apply Particle Splitting to simulations of clump-clump collisions. We investigate the 
influence of clump mass, clump velocity and collision impact parameter on the structures 
formed. We show that such collisions lead to the formation of shocked layers. Networks of 
filaments or spindles, and groups of protostellar discs form within the layers. In all colli- 
sions, fragmentation of the filaments was the common mechanism producing the groups of 
protostars. Low-mass secondary companions to the protostars may form subsequently by 
instabilities in the discs and/or dynamical interaction between the protostars. However, due 
to time-step requirements, we cannot follow the disc evolution for long enough to explore the 
formation of secondaries. 

We show that all the protostars formed have mass accretion rates of ~5 x 10^^ Mq yr~^ 
for the first 10-20 thousand years after their formation. This mechanism shows 10-20% Star 
Formation efficiency. The efficiency increases with increasing clumps mass. Collisions with 
impact parameter b <0.4 and Mach number Ai ^10 give the highest efficiency. We predict 
the existence of filaments with ^5 x 10^ cm~^ in sites of dynamical Star Formation. 
These filaments could be observed in NH3 or CS line radiation, in star forming regions lying 
within 1 kpc. 



Contents 



1 Introduction 



1.1 Stars: the end-product of Star Formation 



1.2 Properties of the ISM 



1.3 Observations of YSOs 



1.4 Th^^ri^ ^ and Simulations of Star Formation 



1.5 Thesis plai 



2 Self- Gravitating Smoothe d Particle Hydrodynamics 



2.1 §^lf-p[r^vitatin^ H y^rodynamic; 



2.2 Fundamentals of SPH 



2.3 Basic SPH equation 
2^ 



2.5 Tree Code Gravity 



2.6 Smoothing length 



2j^^^PI^2ua^ons 



2j^^^itegra^oi^cheme 



2j^^^^d^^^^m^te££ 



2.10 Going through the code 'step bv step 



2.11 Tests 



2.11.1 Colliding Flow! 



2.11.2 Free-Fail collapse 



2.11.3 Stable Isothermal Sphere 



3 Rota ting Clouds wi th m=2 Density Perturbations 



3.1 Initial Conditions 



IX 



3.2 The solution in Eulerian and Laerangian formulations 



^langin^ 



10 



— 1 — 
^ cm 



3j3j2_^Decreasin^^ 



3j3j3_^ncreasing^ 



3.3.4 600.000 particle simulation 



3.3.5 Conclusions 



3.4 Ch anging the number of S PH neighbour 



3.4.1 Isothermal Simulation; 



3.4.2 Simulations with Adiabatic Heating 



3.6 Settled distribution 



4 Particle Splitting 



4.1 Jeans conditio 



1 



4.2 Particle splitting: Concept 



4.3 Particle splitting: Implementatio: 



1 



^3^^^Position^o^tii^fin^£ar^cle£ 



4.3.2 Density profile of the configuration of 13 fine particlei 



^Sj^^^ensit^^tabili^jvi^^articl^gli^ng 



^3^^^moo^in^fing^^o^^^fine_£a£ticte 



^3j^^^elocitie^o^h^fin^£a^cks 



4.3.6 Updating other fluid properties fc numerical parameters 



4.4 Tests 



^^^^^CoUags^^imulation 



4.4.2 Stable Isothermal Sphere 



4.4.3 Rotating Cloud 



4.4.4 Conclusions 



5 Clum p-Clu mp Collisions 



5.1 Initial conditions 



5.2 Repeating the Simulations of Bhattal et al. 



5.2.1 Simulation with Mn=75 MrT^ and 6=0.2 



CONTENTS 





5.2.2 Simulation with Mn=75 and 6=0.4 




5.2.3 Simulation with Mn=75 Mr7^ and 6=0.5 




5.2.4 Conclusions 


5.3 


Low-mass clumo collisions 




5.3.1 Simulation with Mn=10 Mr^. M=5 and 6=0.0 




5.3.2 Simulation with Mn=10 Mr:.. M=10 and 6=0.C 





5.3.3 Simulation with Mn=10 M^. M=15 and 6=0.C 





5.3.4 Effect of chaneine M with constant 6=0.0 






5.3.5 Simulation with Mn=10 Mrr,. M=5 and 6=0.2 




5.3.6 Simulation with /lfn=10 M^. .yU = 10 and 6=0.^ 




5.3.7 Simulation with Mn=10 M^. A^=15 and 6=0.^ 




5.3.8 Effect of chaneine Ai with constant 6=0.2 






5.3.9 Simulation with Mn=10 Mr^.. M=^ and 6=0.4 




5.3.10 Simulation with Mn=10 Mr^.. M=10 and 6=0.4 





5.3.11 Simulation with Mn=10 M,t.. M=15 and 6=0.4 





5.3.12 Effect of chaneine Ai with constant 6=0.4 




5.3.13 Simulations with Mn=10 and 6=0.f 




5.3.14 Simulations with Mn=10 Mrr, and 6=0.^ 






5.4 DiscussiorJ 



6 


Conclusions 








6.1 SPH and Particle SDlittina 






J 

6.2 Cloud-cloud collisions! 



A Jeans criterion of stability 



References 



CONTENTS 



List of Figures 



2.1 


Density and its aradient for an isolated p^irtiglj 


2.2 


1 

Findine SPH neiehbours with the kernel box methodi 


2.3 


Multiole time-steoJ 




2.4 


Tl.e.„H..>,.pHy...p| 




2.5 


SPH tests: Collidine flowsl 




2.6 


SPH tests: Free-fall collaDsJ .... 




2.7 


SPH tests: Stable isothermal sohere 





15 
25 
28 
30 



36 



^J__QQli2nm.d£flSity_ulQli^£-^l£jluMa^^ 43 



3.2 Column density plots for a cloud of 80.000 particles before heating starts and at the end ion = 10 



3.7 Column density plots for a cloud of 600.000 particles before heating starts and at the end ipn = 5 



3.8 Column density plots for the low-resolution isothermal simulations with increasing N„ 55 



3.9 Column density plots for the high-resolution isothermal simulations with increasing N„ 56 



3.10 Column density plots for the low-resolution simulations with adiabatic heating and increasing 



3.11 Column density plots for the high-resolution simulations with adiabatic heating and increasing 



3.12 Density and its gradient for an isolated particle using the kernel of Eqn. H . 65 



3.14 Column density plot for a cloud of 80.000 particles at the end with particles initially taken from a 



4.1 Configuration of 13 fine particles in 2 and 3 dimensions! 76 



4.2 Calculating the density profile of the configuration of 13 fine particles! .... 78 



Xlll 



XIV 



LIST OF FIGURES 



4.4 Gradient of the density profiles of a coarse particle in isolation and of the ensemble of 13 fine particles th 



4.5 Zoom on the density profiles of a coarse particle in isolation and of the ensemble of 13 fine particles that 



4.6 J immediately after the splitting vs. r,- &: the density distribution of the initial settled box (before splitti 



4^^^en^i^^ig^^ii^2^^^h^^^^f^£^gli^n^2^^^^L^/^^^;^^^2ZL 



85 



4^^^en^i^^^roffl^^n^^gn^i^^^^ign^^^^£ng£m^l^£^^^n^£ar^£lg^y^;^^^^ 



86 



4.9 Eyolution of a stable isothermal sphere when nested splitting was applied, without implementation of the 

4.10 The code 'step by step' including particle splitting! 90 



4.11 Radial density profile of a collapse simulation when nested splitting was appliec 



95 



4.12 Radial density profile of a collapse simulation when on-the-fly splitting was applied 96 



4.13 Radial density profile of a collapse simulation when nested splitting was applied with particles initially ts 



^^^^litia^te^^^^teH^sothemia^^ere 



100 



4.15 Eyolution of a stable isothermal sphere when nested splitting was applied . . 103 



4.16 Eyolution of a stable isothermal sphere when on-the-fly splitting was appliec 



104 



4.17 Eyolution of a stable isothermal sphere when nested splitting was applied with particles initially taken fn 
^^^^gg^g^^^li^in^Q^^^Sii^^^^iiffii^Es^isJsS^^SliiSiJ^SiiSi^^^iS^^^^^Si^S^^SS^^S^^^^iJflii^ 



4.19 Nested Splitting for a Cloud of 40.000 Particles: Column density plots of the cloud before heating is a.ppl 

4.20 On-the-fly Splitting for a Cloud of 40.000 Particles: Column density plots after splitting is applied and b' 



4.21 On-the-fly Splitting for a Cloud of 40.000 Particles: Column density plot of the cloud at the enJll4 

4.22 Nested Splitting for a Cloud of 10.000 Particles: Column density plots of the initial sphere and the cloud 



4.23 Nested Splitting for a Cloud of 10.000 Particles: Column density plots of the cloud before heating is appl 
4^^^imida^oi^^^^lou^^^^OjOO^^articlfi^withOT 



4.25 On-the-fly Splitting for a Cloud of 10.000 Particles: Column density plots of the cloud after splitting is a 



4.26 On-the-fly Splitting for a Cloud of 10.000 Particles: Column density plot of the cloud at the endl l21 



5.1 Column density plots for the 6=0.2 collision of two 75 Mm clump 


133 


5.2 Column density plot for the 6=0.2 collision of two 75 Mm clumps showing the network of fllaments produ 


5.3 Column density plot for the 6=0.4 collision of two 75 Mm clumps 




137 
138 


5.4 Column density plot for the 6=0.5 collision of two 75 Mm clumps 




5.5 Column density plots for the 6=0.0 and M=15 collision of two 10 Mm clumps 


145 
146 
149 
150 


5.6 Column density plots for the 6=0.2 and 7W=5 collision of two 10 Mm clumps 


5.7 Column density plots for the 6=0.2 and 7W=10 collision of two 10 Mm clumps 


5.8 Column density plots for the 6=0.2 and A^=15 collision of two 10 Mm clumps 



LIST OF FIGURES 



XV 



5.9 Column densitv plots for the b=OA and A4=5 collision of two 10 Mr^ clumps 153 



5.10 Column densitv plots for the 6=0.4 and M=10 collision of two 10 M^y clumps 154 



5.11 Column densitv plots for the b=OA and 7W=15 collision of two 10 Mrri clumpa 157 



5.12 Column densitv plots for the b=0.6 and A4=5 collision of two 10 clumps 



158 



5.13 Column densitv plots for the 6=0.8 and M=5 collision of two 10 clumpa 161 



LIST OF FIGURES 



List of Tables 



3.1 Summarv of results for the simulations with different on 


51 




3.2 Summarv of results for the low-resolution isothermal simulations with increasing N„ 


54 


3.3 Summarv of results for the hieh-resolution isothermal simulations with increasing 


58 


3.4 Summarv of results for the low-resolution simulations with adiabatic heatine and increasine N„ 


3.5 Summarv of results for the hieh-resolution simulations with adiabatic heatine and increasine N„ 



4.1 Co-ordinates of the 13 fine particles at unit distance awav from each other . . 77 



4.2 Summarv of results for the standard test simulations using particle splitting ffirst Dart) 124 



4.3 Summarv of results for the standard test simulations using particle splitting (second parti 



125 



5.1 Minimum initial numbers of particles & A^™" for different values of Mn and 7^129 



Sj^^^immia^^^imid^on^^n^^ios^n^orten^esid^ 

5.3 Dependence of different Quantities and phenomena on the increasing values of b and M for simula 



xvii 



xviii LIST OF TABLES 



Chapter 1 



Introduction 



Star Formation is a field relying on theoretical work since it is difficult to obtain observa- 
tionally a picture of the way stars form. Numerical simulation of Star Formation has lately 
become a very popular theoretical tool due to the rapid increase in computing power. In this 
thesis, we deal with simulations of low-mass star formation and in particular, we establish a 
method for increasing the resolution and dynamic range of simulations. In this chapter, we 
briefly review the main theoretical and observational constraints on Star Formation, so that 
we can put our work in context with the other research in this field. 

1.1 Stars: the end-product of Star Formation 

Stars are a fundamental constituent of the Universe. Stars are hot massive dense gas spheres 
emitting radiation produced in their centres from nuclear reactions and/or gravitational con- 
traction. Stars are very well studied. Analysing their spectra provides information not only 
on their outer layers that can be seen emitting, but also on their interiors. In particular, we 
can infer their chemical composition and from known nuclear reaction cycles we can produce 
models for their structure and time evolution. Core hydrogen burning stars of different age, 
mass and chemical composition define a locus on the surface-temperature vs. luminosity 
diagram: the Main Sequence of the Hertzsprung-Russell diagram. A star reaches the Main 
Sequence, following a period of gravitational contraction, as soon as the conditions at its cen- 
tre become hot and dense enough to start burning hydrogen. The gravitational contraction 
phase of a forming star is qualitatively understood and there are theoretical models predict- 
ing its evolution towards the main-sequence (pre-main-sequence tracks). What is relatively 
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unknown, and constitutes one of the fundamental questions for Star Formation theory, is the 
mechanism that brings star-forming gas together in the first place. At this point, a review of 
the properties of the interstellar medium, which provides the ingredients for Star Formation, 
can help us put this question into perspective. 



1.2 Properties of the ISM 

The Interstellar Medium (ISM) consists of gas in all states (atoms, molecules, ions) and dust 
grains. The dust component accounts only for 1-2% of the mass of the interstellar medium. 
The gas consists mainly of hydrogen and helium. A small percentage of the gas mass is in 
the form of heavier elements. 

The ISM is far from being uniform, as it consists of dense clouds of gas {n > 100 cm^'^). 
In these clouds, the gas is mainly molecular and cold (T < 100 K). Apart form II2, the 
clouds contain quite a few other molecules, such as CO, OH, CH, II2O, CS, NII3 as well 
as more massive and more complex carbon chains ( |van Dishoeck &: Hogerheijde 1999 1. Such 



clouds are called Giant Molecular Clouds (GMCs) because of their size (R ~ 10-20 pc and 
M ~ lO'^-lO'^ Mq). These clouds are immersed in a warm (T > 10^ K) and rarefied {n > 1 
cm~^) interclump medium, partly atomic and partly ionised. 

The structure of the ISM is evolving with time. Parts of a GMC in proximity to massive 
stars can be ionised by expanding nebulae (HII regions, supernova remnants) and become 
part of the interclump medium, or be squashed and forced to form more stars. GMCs 
form in galactic spiral arms possibly through condensation of HI clouds and /or collisional 
agglomeration of lower mass clouds. GMCs live for about 10^ yr. 

GMCs are self-gravitating and they are supported by supersonic MHD turbulence. Some 
of the turbulent motions are generated by stars within these clouds. In particular, stellar 
winds and outflows, HII regions and supernova remnants are mechanisms for exciting tur- 
bulent motion. There is also evidence for magnetic fields threading the clouds (observed via 
polarimetry and Zeeman splitting measurements, e.g. Ward-Thompson et al. (|2000|) ). 

Molecular Clouds are the birthplaces of all known young stars. They provide the initial 
conditions for Star Formation. "They host Young Stellar Objects (YSOs) in a wide range 
of evolutionary states; from Class protostars some 10-2 Myr old, der iving most of their 
luminosity from gravitational infall, to T-Tauri stars a few Myr old, deriving their luminosity 
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from quasi-static contraction. They also host stars in a wide range of spatial groupings; from 
isolated single stars as in Taurus, having no known neighbours within a few pc, to rich star 
clusters as in Orion, having a few thousand stars within a few pc. The masses of the stars 
in GMCs range from 0.1-30 M©, nearly the whole range of known stellar masses. Indeed the 
Initial Mass Function (IMF), or the distribution of stellar masses at birth, is indistinguishable 
between stars in molecular clouds and field stars." ( [Myers 1999 ) 

GMCs have structure of their own. Dense clumps can be observed within GMCs, typ- 
ically with densities above 10'^ cm~'^, via molecular line observations using tracers such as 
CO and its isotopomers. These clumps have masses between 10 and 100 Mq. In sites 
of high-mass star formation, these clumps tend to be associated with groups and clusters 
of young stars (ILada et al. 199l|l . Within these clumps, dense cores can be observed using 
molecules such as CS and NH3 which trace densities of 10^ cm~'^ and above ( Myers et al. 1983) 
Myers &: Benson 1983| [Myers 1983"! [Benson fc Myers 1983D . In sites of low -mass star forma- 



tion, such cores (with masses of ~1 Mq) are directly associated with the formation of single, 
or at the most binary, stars. It is, therefore, of great importance to understand how these 
cores and clumps can evolve to produce young stars. Before we list some of the available 
models for Star Formation, let us summarise the properties of young stars. 



1.3 Observations of YSOs 

Some of the cores observed in NH3 are associated with IR sources which are identified as 
protostars. There are also starless cores or pre-stellar cores ( [Ward-Thompson et al. 1994 ), 



which are believed to be the precursors of protostars. They are more extended than cores 
containing IR sources. Some show double-peaked velocity profiles with a stronger blue-shifted 
peak suggesting that they are collapsing ( [Ward-Thompson et al. 1996 1 . Pre-stellar cores are 



believed to be supported by magnetic fields. The magnetic support is removed by ambipolar 
diffusion (jShu et al. 1987|1 or by reconnection of magnetic field lines ( Lubow &: Pringle 1996 ) . 
Without this support the pre-stellar cores start collapsing. 

The cores associated with IR sources also show signatures of collapse. The YSOs within 
these cores are not all at the same evolutionary stage. The evolution of YSOs can be divided 
into two major phases: the embedded phase and the revealed phase ()Lada 1999|) . There are 
two different classes of objects in each phase. 
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In the embedded phase the protostars are collapsing fast. Their luminosity is produced by 
gravitational contraction and accretion. They are characterised either as Class or as Class 
I objects depending on their Spectral Energy Distribution (SED). Class objects have SEDs 
that can be fitted with single temperature black body functions. The Class SEDs peak at 
sub-mm wavelengths and the objects are not observed below 20 fim. Class protostars are 
highly obscured (Ay ~ 1000) and of extremely low temperature (20-30 K). Class SEDs 
are attributable to emission from dust grains in the infalling envelope. The dust absorbs 
all radiation coming from the protostar and then re-emits it in the sub-mm spectral range 
()Andre et al. 1993|l . The Class phase lasts for about 1-5 x 10^ yr. During this period, Class 
protostars accrete matter with a mass accretion rate of > 10^^ Mq yr~^. These objects 
are associated with very powerful outflows. 

The SED of a Class I object is broader than a single temperature blackbody function. 
Class I SEDs peak at sub-mm wavelengths but they also show a characteristic excess of 
infrared emission. This IR emission is associated with large amounts of circumstellar dust. 
Class I sources are deeply embedded and there is no emission at optical wavelengths. However, 
there is NIR emission (e.g. 2.2 /im) and a significant fraction of this NIR emission is from 
scattered light. Class I objects are also associated with outflows, but these are not so powerful 
as those of Class objects. The lifetime of a Class I object is 1-5 x 10^ yr. During this time, 
the mass accretion rate is < 10^^ yr^^. There is some evidence that Class I objects have 
come through a Class phase. 

In the revealed phase the protostars have evolved into pre-main sequence stars. The 
infalling envelope of gas and dust has been removed. The luminosity in this phase is pro- 
duced by Kelvin-Helmholtz contraction and deuterium burning. Again, there are two classes 
of objects: Class II and Class III. Classification of pre-main-sequence stars into these two 
classes is based on their SEDs. The Class II SEDs are broad like those of Class I objects 
but they peak at IR wavelengths. For wavelengths longer than 2.2 ^um they decrease in a 
power-law fashion. The IR excess of Class II objects is smaller than that of Class I objects. 
Again, it indicates the presence of circumstellar dust. These YSOs can be observed in optical 
wavelengths, and therefore are better known than the protostars of the embedded phase. 
They show Hq, emission which is associated with accretion onto the protostar. Optical iden- 
tification of accretion discs silhouetted against a bright background has been possible with 
HST ()0'Dell fc Wen 1994|) . Some of these objects have weak outflows associated with their 
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accretion discs. Class II objects have variable emission lines; they are also known as Classical 
T-Tauri Stars (CTTS). From pre-main-sequence tracks, their age is estimated to be 1-4 x 10^ 
yr. Their accretion rates are estimated to be ~ 10^® Mq yr^"*^. 

During the infall of circumstellar gas onto a protostar, its angular momentum increases. 
Above a critical value, the excess in angular momentum, produced by more infalling mass, 
has to be removed. The discs provide a natural mechanism for angular momentum removal 
( Lynden-Bell &: Pringle 1974 ) . Mass is accumulated on the plane perpendicular to the angu- 
lar momentum vector and angular momentum is transfered to the mass at the outer edge of 
this disc. 

Class III objects are visible at both optical and IR wavelengths. Their SEDs can be 
fitted with single temperature black body functions. At wavelengths longer than 2.2 ^um 
their emission decreases more steeply than that of Class II objects. They show no excess IR 
emission. Their SEDs are thought to be arising directly from their photospheres. However, 
there is still considerable amount of dust around these objects (remnant of the infalling enve- 
lope of previous protostellar phases) which produces a substantial extinction and reddening. 
These objects could not be distinguished from main-sequence stars, if they were not emitting 
hydrogen lines and X-rays. The Hq, emission is not so strong as in the Class II stage, but 
it indicates traces of a disc. Class III objects are also known as Weak-line T-Tauri Stars 
(WTTS). Their positions on the H-R diagram can be directly compared to predictions of 
pre-main-sequence tracks (they are positioned on top and to the right of the main-sequence). 
From these tracks we obtain their ages of 10^ - 10^ yr. 

The discs around Classical T-Tauri stars are rather extended (100-1000 AU) but not very 
massive (0.001-0.01 Mq). The central protostar, which has accreted most of its mass by this 
time, has mass 0.5-1.5 Mq (jBeckwith 19991) . At the end of their evolution (Class III stage), 
pre-main-sequence stars have their discs stripped of their gaseous component. What little 
gas remains, together with the dust, is believed to form planetary systems ()Ruden 1999|) . 

T-Tauri stars were first observed in the Taurus molecular cloud. Taurus is a low-mass 
Star Formation Region (SFR) on the northern sky at a distance of ~140 pc. Taurus is the 
best example of distributed star formation, where forming stars are not part of dense groups 
or clusters (|Gomez et al. 19931 IGladwin et al. 1999|) . Young stars in Taurus have masses 
between 0.5 and 1.0 Mq and are in small clusters of <10 members. 

The Orion molecular cloud (also on the northern sky, at a distance of ~440 pc) is the 
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most well-studied example of clustered star formation. It has been shown that in the center 
of Orion (close to the Trapezium OB stars) there are more than a thousand stars within one 
or two pc ()Hillenbrand 19971 IHillenbrand &: Hartmann 1998|) . In Orion B, almost all stars 
(~700) have formed in just 4 clusters (jLada et al. 199l|) . Discs around young stars in Orion 
(proplyds) have smaller radii than those in Taurus possibly due to photo evaporation of gas 
by the radiation field of the nearby massive stars ()Hollenbach et al. 1994|) . This illustrates 
the strong effect that the molecular cloud environment has on Star Formation (jLada 1999|) . 

The more massive of the YSOs in Orion (M >2 Mq) are Herbig Ae-Be stars and are the 
precursors of main-sequence stars of type A or earlier. These YSOs are not so well studied 
as T-Tauri stars because they are fewer of them and they are further away. Despite having 
luminosities large enough to drive stellar winds, these YSOs are highly obscured by their 
infalling envelopes. Like T-Tauri stars, they are associated with violent phenomena like jets 
and outflows (|Churchwell 1999|) . 

It has been shown that the IMF, (j){M), can be fitted with 



(p{M) dM ~ 50 exp 



log?o 



M 



O.2M0 ^ M ^ 5OM0 



O.OIMq 

(|Miller fc Scalo 1979jl . This IMF extends the previous Salpeter H1955jl IMF 

0(M) dM ~ KM-^-^^dM, OAMq < M ^ IOM0 

to smaller and larger masses. In the low-mass range the Miller & Scalo IMF is flatter than the 
-2.35 power law, whereas in the high mass range it is steeper than the -2.35 power law. This 
IMF suggests that low-mass stars are greater in number than high-mass stars, and account 
for most of the mass. 

The observed IMF puts a strong constraint to Star Formation theories, as these theories 
must predict a similar mass distribution. Another statistical constraint is set by the observed 
spatial distribution of stars. 

It is well known that almost 50% of all solar-type main-sequence stars are part of a binary 
or higher multiple system ( Duquennoy Hi Mayor 199T I . Duquennoy & Mayor have measured 



the distributions of periods, eccentricities and mass ratios for field binary systems. They 
have shown that the distribution of periods can be fitted with a Gaussian that peaks at ^^200 
years, corresponding to a separation ~30 AU. Pre-main-sequence binaries are observed in 
the same range of separations as field star binaries (from a few to 10^ AU, e.g. Ghez et al. 
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H1997|) ). Their distribution of separations is similar to that of the field stars (e.g. Prosser et 
al. H1994|) ^. This indicates that stars often form in binaries. 

The binary frequency may depend on the environment of the SFR (jBrandner &: Kohler 1998|l . 
It appears that closer binaries are formed in regions near massive stars. However, there is 
evidence that universal processes determine the multiplicity of young stars, since the surface 
density of companions for pre-main-sequence stars can be fitted with a single power law, 
N{e) oc 6i-2.o±o.i fQj. 0.0001° ^ 6* ^ 0.01°, in many different regions, e.g. Chamaeleon I, Tau- 
rus, Orion, p Ophiuchus, Lupus, Vela (Larson (|1995|) : Simon H1997|) : Nakajima et al. H1998|) : 
for a review see Gladwin et al. H1999|) ^. 

We have seen that protostars, which are very young but not very bright sources, are 
heavily obscured by circumstellar material for a significant period of their formation. It 
is, therefore, rather difficult to obtain detailed observational information for the physical 
mechanisms that convert the ISM into stars. Theoretical work therefore has an important 
role to play in Star Formation research. Modelling the complex gas dynamics in the ISM, 
including several orders of magnitude change in density and linear scale, is achieved by means 
of numerical simulations. In this thesis, we model cloud-cloud collisions. In the next section, 
we review the most important theories of Star Formation and previous cloud-cloud collision 
simulations. 

1.4 Theories and Simulations of Star Formation 

"Every piece, or part, of the whole of nature is always merely an approximation to the 
complete truth, or the complete truth so far as we know it. In fact, everything we know is 
only some kind of approximation, because we know that we do not know all the laws as yet. 
Therefore, things must be learned only to be unlearned again or, more likely, to be corrected." 
( Feynman 1995 ). 

Star Formation theories should predict the formation of groups of protostars with the 
above distributions of masses and separations. One of the oldest mechanisms proposed for 
the production of binary stars was fission, the splitting of a single object into two due to 
excess rotation. However, both finite difference and SPH simulations of Durisen et al. (|1986|) 
have shown that spiral arms remove angular momentum efficiently from the rotating object, 
and fission does not occur. 
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Capture predicts that an initially unbound system of two or more protostars will be- 
come bound when it loses energy in a close encounter. Star-disc capture ()Hall et al. 19961 
IBoffin et ffl/."l9 981 is a mechanism sufficient to reproduce the binary frequency in very small 
dense stellar clusters, but not in larger looser clusters or the field ( Clarke Hi Pringle 1991 1 . 

The well known model of low-mass star formation of Shu, Adams & Lizano H1987() involves 
the collapse of an isothermal sphere producing a single protostar. The isothermal sphere loses 
its magnetic support via ambipolar diffusion, where neutrals slowly decouple from the ions 
and the field to produce a p oc density profile. Collapse of a sphere with such a profile 
gives a constant accretion rate. However, Whitworth et al. H1996() have argued that such 
a profile is unlikely to arise in nature. In addition. Class protostars are observed to be 



undergoing collapse with a less centrally condensed proffie ( Ward-Thompson et al. 1994 ). 

Fragmentation is a process where separate parts of a cloud, clump or core become grav- 
itationally unstable and contract. Groups of protostars are usually produced by this pro- 
cess. Fragmentation requires a mechanism to induce gravitational instabilities. Several 
numerical simulations of fragmentation have been conducted involving clouds of different 
geometries or applying different mechanisms to produce the instabilities (for a review on 
fragmentation simulations see Bonnell (|1999j) ). Simulations of rotating spherical clouds 
()Boss fc Bodenheimer 19791) or cylindrical clouds (jBonnell fc Bastien 1992} IBonnell et al. 1991]) 
produce binaries and/or bars that break up into lumps. In disc fragmentation simulations, 
spiral arm rotational instabilities in the rotating disc produce companions to the central 
objects (|Bonne]] 1 9941 IBonnell Bate 19941) . 

Mechanisms of a dynamic nature that can induce gravitational instabilities and binary 
formation via fragmentation, include the creation of shock compressed layers triggered by the 
interaction of stellar winds with the ambient gas, or by cloud-cloud collisions. 

Early simulations of cloud-cloud collisions produced coalescence of clouds (|Stone 1970a| 
IStone 1970bl IGilden 19841 ILattanzio et al. 1985|) . Fragmentation of a shocked layer produced 
at the interface of the collision was found in simulations with equations of state that included 
cooling of the shocked gas ( Nagasawa Miyama 19871 IHahe »hta 19921) . 

Subsequent SPH simulations of cloud-cloud collisions ( Chapman et al. 1992} Pongracic et al. 19921 
ri'urner et al. 19951 IWhitworth et al. 19951 IBhattal et aTm )^) have been able to resolve sev- 
eral orders of magnitude change in density and linear scale. These simulations produced 
large rotationally supported protostellar discs (200-4000AU) of high mass (5-60 M©). The 
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simulations of Chapman et al. (|1992|) involved higher mass clouds and produced filaments 
that fragmented, producing several discs. Rotational instabilities in the discs produced sec- 
ondary companions to the protostars (jTurner et al. 19951 IWhitworth et al. 1995|) . Increasing 
the impact parameter of the collision ( Pongracic et al. 1992| iBhattal et al. 19 98 ) assisted disc 
fragmentation as the angular momentum was increased. However, such simulations did not 
always obey the Jeans condition (jTruelove et al. 19971 IBate Burkert 1997|) and as a result 
they are suspect, i.e. the protostars may have formed numerically. 

In this thesis, we apply a new method (particle splitting) that enables the Jeans condi- 
tion to be obeyed throughout our simulations with minimum computational cost. With this 
method, the resolution of the simulations can be increased at certain times and places, cre- 
ating simulations with increasingly fine resolution, nested inside an initial coarse simulation. 
We apply particle splitting to numerical simulations of cloud-cloud collisions. 



1.5 Thesis plan 

In chapter [21 we describe in detail the numerical code we use (Smoothed Particle Hydrody- 
namics with Tree-Code-Gravity). Most of the features of this code have been discussed in 
previous theses ( [Chapman 19921 ri'urnerl993l IBhattal 19961 IWatkins 1996|1 and in Turner et 
al. p995|) . but they are presented here in detail in order to make the thesis self-contained. 
At the end of the chapter, we apply our code to standard tests for a self-gravitating hydro- 
dynamical code, to demonstrate that it is appropriate for the problems we study. 

In chapter |31 we present simulations of rotating clouds with an m=2 density perturba- 
tion. With these simulations, we try to reproduce previous results (|Bate fc Burkert 19971 
ITruelove et al. 19971 ITruelove et al. 19981 IKlein et al. 1999j) . and we investigate the depen- 
dence of the code on numerical parameters such as the number of particles, the number of 
neighbours, the choice of the interpolating function, the initial spatial distribution of the 
particles. We consider this chapter to be an examination of the efficiency of the code when 
applied to more realistic problems. 

In chapter|3 we describe particle splitting in detail. The method is then tested thoroughly. 
From the tests, it is shown that the noise introduced to the simulations when the method is 
applied, is not significant. We also apply particle splitting to simulations of rotating clouds 
with m=2 density perturbations: results of finite difference simulations are reproduced with 
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great computational efficiency. 

In chapter |SJ we apply particle splitting to simulations of cloud-cloud collisions. We 
repeat the simulations of Bhattal et al. (|1998|) but now satisfying the Jeans condition. These 
simulations involve collisions between intermediate-mass clouds (75 M©). By comparing our 
results to those of Bhattal et al, we can estimate the efficiency of the new method and 
quantify its benefits. We then present simulations of low-mass clump collisions (10 Mq) that 
have not been studied before. We investigate the influence of cloud mass, collision impact 
parameter and relative velocity on the filamentary structures formed in the shocked layers, 
on the protostellar discs formed within the ffiaments and on the Star Formation efficiency. 
Finally, we compare the properties of the stellar systems produced in our simulations and 
those obtained from observations of YSOs. 

In chapterlHJ we summarise our main conclusions emphasising the computational efficiency 
achieved with particle splitting. Finally, we make suggestions for improvements to our code 
as well as for further applications of particle splitting. 

Appendix A gives a derivation of the Jeans criterion for fragmentation. 



Chapter 2 

Self-Gravitating Smoothed Particle 
Hydrodynamics 



2.1 Self-gravitating Hydrodynamics 

The gas in the interstellar medium is highly compressible. In star formation the self-gravity 
of the gas is also important. In order to describe the evolution of a self-gravitating, inviscid, 
compressible, non-magnetic fluid, we need to solve a system of four equations - the continuity 
equation, Euler's equation, energy equation and equation of state ()Landau &: Lifshitz 19661 
IShu 1992]) - with four unknowns, namely the velocity v, pressure P, specific internal energy 
u, and density p, at each position r in the fluid. It is implicit that these four quantities are 
also functions of time. The four equations read as follows: 

• Continuity equation 



where we have used the Lagrangian time-derivative d/dt = d / dt + v ■ V and the 
fact that dp / dt = — V • {pv). The continuity equation expresses the conservation of 
mass. 

• Euler's equation 
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where agrav is the self-gravitational acceleration, given by 



J all space I '■I 

and a^isc is the artificial viscous acceleration (see discussion in ^2.4|) . Euler's equation 
expresses the conservation of momentum. 

• Energy equation and equation of state 

In general, the equation for the rate of change of the specific internal energy reads 



P(r)^ = -P(r)V-^(r) + (E - A), 

where E and A are the radiative heating and cooling rates per unit volume respectively. 
This equation expresses the conservation of energy. The pressure is then given by the 
ideal gas equation of state 

P = (7 - l)pu, 
where 7 is the ratio of specific heats. 



However, it is well established (jTohline 1982|) that prestellar gas at low densities {p <C 
po ~ 10^^'* g cm~'^) is approximately isothermal at T ~ lOK (this is mainly due to the strong 
temperature sensitivity of A). At high densities (p ^ 10^^^ g cm^^), the gas is approximately 
adiabatic. Therefore, we can reduce the system of equations we have to solve, by using a 
barotropic equation of state instead of the two equations involving the specific internal energy. 
The barotropic equation of state which we use reads 



P(r) 



1 + 



Po 



4/3' 



11/2 



(2.4) 



where cq is the isothermal sound speed of the gas at low densities and po is the density above 
which the gas becomes adiabatic. Eor p <C po Eqn. 12.41 gives P ~ pc^, while for p ^ po, 
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At high densities, the gas is approximated with a polytrope, having an index of 3/2 ^, 
only for a few orders of magnitude. For further evolution toward stellar densities detailed 
radiation transport calculations are necessary. To-date, no simulation in the literature has 
advanced that far. Bate (^98) has presented results approaching stellar densities, but with 
an approximate barotropic (piecewise polytropic) equation of state. 

With our code, we do not attempt to model interactions between the gas and interstellar 
magnetic fields that thread molecular clouds. This would require a larger set of governing 
equations (e.g. see the MHD equations in Vazquez-Semadeni, Canto & Lizano (|1998|1 ). but 
more importantly it would require a two-, or possibly a three-, component fluid. Such a 
fluid has never been modelled with an SPH code and it is unknown what the resolution 
requirements for each fluid component would be. Exploring this is beyond the scope of the 
present thesis; and in any case, it is not clear that the magnetic fields observed in molecular 
clouds are sufficiently strong to greatly influence the dynamics of the cloud-cloud collisions 
studied in this thesis. 

To solve the above system of equations, we have used the numerical method described 
in Turner et al. (|1995|) with some later improvements. The method consists of two numer- 
ical techniques: SPH, which provides values for the hydrodynamical properties of the fluid 
(3231), and Tree-Code-Gravity (TCG), which provides values for the gravitational accelera- 
tions f ^2.5() . In this chapter, we shall describe in detail the method as well as the integration 
scheme with which we follow the evolution of the fluid in time f ^2.8() . To summarise the chap- 
ter, we will give a schematic picture of a complete cycle of the integration scheme f' j2.10j) . 
Finally, we will briefly describe the performance of the numerical method when applied to 
some standard tests for a self-gravitating hydrodynamical code f ^2. 11(1 . The novel modifi- 
cations to this numerical method associated with particle splitting are developed in detail 
in chapter [l] Additional numerical techniques used for the purposes of special applications 
of this numerical method, are described in the relevant chapters, e.g. radiative cooling of 
shocked layers in chapter El sink particles in chapter^ 



^The two rotational degrees of freedom that bring Tohline's 1198211 polytropic index to 5/2 are frozen out 
while the gas is below ~400K ijWinkler fc Newman 198?)|l . Above this temperature, the rotational degrees of 
freedom are excited, and the 5/2 index switches on for one or two orders of magnitude increase in density, 
before dissociation happens and the gas finally becomes atomic with an index of 3/2. 
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2.2 Fundamentals of SPH 



Smoothed Particle Hydrodynamics (SPH) ( [Gingold Monaghan 1977| Lucy 1977 1, is a La- 



grangian numerical method that assumes no symmetries or imposed grids. It is, therefore, 
very efficient in describing problems which involve complex 3-dimensional geometries. SPH 
represents the fluid by discrete but extended/smoothed particles (i.e. Lagrangian sam- 
ple points). The particles are overlapping, so that all the physical quantities involved can 
be treated as continuous functions both in time and space. To implement this, a smooth- 
ing function (kernel) with compact support is used. This smoothing function describes the 
strength and extent of a particle's influence. Bhattal H1996() has shown that the M4-kernel, 
that has been frequently used in SPH ( Monaghan fc Lattanzio 1985 1 , gives good results. The 



3-D M4- kernel is a polynomial 



Wma{s) = - 



1 - 3sV2 + 3sV4, ^ s ^ 1; 

(2-s)3/4, l^s^2; (2.5) 

0, s ^ 2, 



which smoothes the mass of a particle out to two smoothing lengths. In SPH, the value of 
any quantity A at a position r is evaluated using 



A{r) = Y.^^^K'w(^-^), (2.6) 



where is the position, rrii the mass and hi is the smoothing length of particle i ( Monaghan &: Lattanzio 1985} 



Monaghan 1988[ Monaghan 1992 ). Ai and pi are the values of A and p at rj. Note the nor- 



malising term which is introduced when we use Eqn. 12.51 for the kernel and substitute 
s = \r\/h. In Eqn. 12.61 the summation is finite due to the fact that the kernel function has 
compact support. This implies that contributions from only a few close neighbouring particles 
are taken into account. The fact that the summation does not need to be over all particles 
greatly reduces the computational cost of all SPH calculations. The benefit of this fact must 
be balanced against the need to obtain accurate results. With h being large enough we can 
reduce sampling errors. The value of h is, therefore, of great importance. It is discussed in 
detail later in this chapter. 
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Figure 2.1. Left : Radial density profile of an isolated particle of unit mass and smoothing length h. 
Right: The gradient of the radial density profile of an isolated particle of unit mass and smoothing 
length h. 

For example, the density at position r, is given by 



p{r) = '^niih- 



r - Ti 



hi 



(2.7) 



The density profile of an isolated particle of unit mass and smoothing length h, calculated 
with Eqn. \2.7\ is shown on the left panel of Fig. 12.11 

The gradient of any quantity yl at a position r is evaluated using 



V^(r) = ^mi—h-^W' 

Pi 



> r -n \ r - n 



r - Ti 



(2i 



where W' (s) = dW/ds ( |Monaghan Lattanzio 1985| [Monaghan 1988} |Monaghan"T992 1 . 

Note that in such an expression, a second order truncation error^, 0(/i^), is introduced 
()Moore 1995|) . For the M4- kernel, the gradient is given by 



3s - 9sV4, ^ s ^ 1; 
Wuiis) = I 3(2 - 5)2/4, 1 ^ s ^ 2; 



(2.9) 



Tliese errors axe introduced because we substitute tlie integral expressions witli sums jMonaghan "19921 1. 
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The gradient of the density is therefore given by 

Vp(r) = J2m,h-'W' ^ (2.10) 

i 

The gradient of the density profile of an isolated particle of unit mass and smoothing length 
h, calculated with Eqn. 12.101 is shown on the right panel of Fig. 12.11 



2.3 Basic SPH equations 



Following Eqn. 12.81 the first two of the set of three hydrodynamical equations that we need 
to solve (namely Eqns. 12.11 &: [T^ should read as follows: 

^-''«E^-'».-'(^)^ (-) 

and 



—r— = mihi^—W — ]- r + 8igrav{r) + Sivisc{r), (2.12) 

dt p{r) ^ ^ Pi \ hi ) \\:- r^l 

where s^grav is given by Eqn. 12.31 and a^jj^c is discussed in ^2.41 

However, in such a case, we would have to use the value of /o(r), numerically estimated 
using Eqn. 12.61 which would introduce an error greater than the second order truncation 
error ()Moore 19951). Instead, we use a formulation that includes density in the gradient of 
the function. This formulation derives from the identities 

pSJA = V(p^) - 

and 

VA ^,A^ A^ 

= V(-) + -Vp. 

p P 

Substituting for v and P respectively, we obtain 



pV ■ V = V ■ {pv) — V ■ Vp, 



(2.13) 



and 



VP ^,P, P^ 

= V - + -Vp. 

P P P^ 



(2.14) 
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dp(r) 



mi/i- (■i;(r) - Vi) ■ W 



(2.15) 



and 



''^'^ - E-A-^(| + ^)^'(V^)^+-.^"^«+Wr). (2.16) 



dt 



As mentioned at the beginning of the previous section, SPH fohows the evolution of the 
hydrodynamical properties of a fluid represented by a system of particles - sample points. 
These particles are allowed to move with the fluid, as they trace elements of constant mass. 
The particle motion will be discussed further in ^2.61 A consequence of this fact is that 
we only need to estimate the fluid hydrodynamical properties at the particle positions, j. 
Therefore, Eqns. 12.151 &: 1^.161 reduce to 



and 



where Vij = Vj — Vi, Vij = Vj — and hij = 0.5(/ii + hj). The latter is used in order to utilise 
the symmetrical forms that Eqns. 12.171 &: [2.181 have obtained with the use of Eqns. I2.1HI fc 
12.141 This way we make sure that every interaction between pairs of particles is symmetric 
and hence we ensure that linear and angular momentum are conserved. 

To calculate density we use Eqn. 12.71 instead of the continuity Eqn. 12.171 This reads as 



This formulation conserves mass very accurately because of the fact that the kernel function 
has compact support and it is appropriately normalised. It is therefore prefered to Eqn. 12.171 
It is also simpler and faster to calculate. In our integration scheme f ^2.8|l . we use Eqn. 12.181 
to calculate the acceleration. 
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The equation for the rate of change of the specific internal energy can also be formulated 
in a symmetric form (Monaghan 19921, 



dt 



1 



rriih, 




Vi, ■ W 



<-i] I 



hi 



H ruih 

Pi V 



-3 ( r^ ^^\ y I I'-'^ji 



pi 



hij 



However, here we are using a barotropic equation of state (Eqn. 12.41) to obtain the pressure. 
In our formulation, the barotropic equation of state reads as follows: 



12 
Pj 



1 + 



El] 

PoJ 



4/3- 



1/2 



(2.20) 



2.4 Artificial viscosity 



In regions where particle streams collide interpenetration may occur. This is undesirable since 
colliding fluid elements should retain their relative positions; the colliding streams should be 
decelerated by shocks. 

In order to make sure that particle interpenetration is inhibited and that we obtain well- 
defined (but not necessarily well-resolved) shocks, we have to include an 'artificial viscosity'. 
We use the artificial viscosity described in Monaghan ()1992|) . This includes a linear bulk 
viscosity component that prevents interpenetration as well as a Von Neumann-Richtmeyer 
type viscosity component. The artificial viscous acceleration that acts on a particle j is given 
by 



mih-;'u,,w 



Z^'— ij U V 



(2.21) 



where 



Pi] 



0, 



, {vij ■ Tij) < 0; 

{Vij ■ Tij) > 0, 



(2.22) 



and 



{Vij • Tij'jhij 



. .12 



+ 0.01/i?, 



(2.23) 



while pij = 0.5(/3j + pj) and Cij = 0.5(ci + Cj) (the average sound speed). The Vij ■ Yij < 
condition makes sure that only particles that are approaching particle j will contribute to its 
artificial viscous acceleration. The above formula is symmetric in order to make sure that 
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particle j will in turn exert an equal and opposite artificial viscous acceleration on any of the 
neighbouring particles, i, that are approaching. 

Hij has dimensions of velocity squared over density. The a-term gives an artificial viscous 
acceleration similar to the hydrodynamical acceleration if sound speed squared is replaced by 
the product of sound speed times the relative velocity. It can treat a bulk colliding flow, but 
in high Mach number shocks the /J-term is needed, as in such cases, the relative velocity of the 
colliding flows is large compared to the local sound speed, a and (3 are tunable parameters 
and should take appropriately large values to prevent interpenetration. If a = /? = then 
the artificial viscous acceleration is zero, particles penetrate each other and shocks cannot 
be modelled. If a and (3 are very large then the shock is very broad and it ends up not 
very well resolved; even mild sound waves are rapidly damped. We have adopted the value 
a = (3 = 1. It has been shown that these values give good results f ^2. 11.11 also Patsis (1999 
private communication)). 

Watkins et al. (|1996)1 have shown that one can use the Navier-Stokes equations for viscous 
compressible flows to derive a viscous acceleration similar to the a term, but the [3 term is 
still needed to model accurately high Mach number shocks. The Navier-Stokes viscosity was 
derived to give a more realistic treatment of shear viscosity (e.g. in problems involving the 
evolution of protostellar discs). The formulation of Watkins et al. H1996|) is derived from the 
standard SPH cross product expressions (Monaghan 1992 1 . 



The artificial viscous acceleration of Eqn. 12.211 should be added to the hydrodynamical 
acceleration of Eqn. I2.1(S| as we shall see in ^2.71 

2.5 Tree Code Gravity 

We now hctv6 to define tlie (igrav^j 

term of Eqn. 12.181 In the case of point masses, the 
gravitational acceleration on each particle j should be calculated as 

^grav,j — ^ ^ ^i'l Tg; (2.24) 

where we have used units so that G =\. However, for a large number of particles, A'^, such a 
formulation becomes very expensive, as the computational cost scales as 0{N'^). Instead, we 
implement Tree-Code-Gravity (TCG) (|ljarnes fc Hut 19861 [Hernquist Katz 1989| ), which 
scales as 0{NlogN). 
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With this method, a tree is constructed containing spatial information on individual 
particles and the centres of mass of groups of particles. This way, for distant interactions we 
can substitute individual particles with groups of particles. At every integration we have to 
construct the tree, walk up the tree to calculate the centres of mass and walk down the tree 
to calculate the gravitational acceleration. 

For the construction of the tree, we use the whole computational domain as the rootcell, 
the top level of the tree. We then divide the rootcell into 8 subcells. These subcells define 
the first level down the tree. Each of them is subsequently divided into another 8 subcells, 
i.e. the next level down in the tree, etc. A cell at any level is not divided further only when 
it contains either a single particle or no particle at all. 

If a cell contains more than one particle, its centre of mass is calculated using all its 
subcells in all lower levels. This way, for the rootcell we calculate the centre of mass of the 
whole computational domain. For each cell, we save the following information: cell centre, 
linear dimensions, pointers to its subcells, total mass and centre of mass, pointers to particles. 

We calculate the gravitational acceleration at particle j using the centre of mass of a 
cell unless the cell is so close that we must use its subcells instead. We apply a geometrical 
criterion in order to decide whether we shall use a cell or its subcells. Specifically, a cell is 
used if 

where I is the linear size of the cell under consideration, D is the distance between particle j 
and this cell and 6 is the maximum opening angle, an accuracy parameter. If the criterion 
is not satisfied then the cell is 'opened' and its subcells are examined. If a subcell contains 
a single particle then a particle-particle interaction is calculated. The value of 9 should be 
sufficiently small for close interactions to be calculated as particle-particle. Salmon, Warren 
& Winckelmans (I1994j) examined different values for 9 and they determined that accurate 
calculations of the gravitational acceleration were obtained for 9 < 0.577. We have chosen the 
value of = 0.576 again driven by the need to balance accuracy against speed of calculation. 

Therefore, the gravitational acceleration of particle j is the sum of contributions from 
other particles and cells. For a particle-particle interaction with particle i we use 

^ij 

I ^ij I 



(2.25) 
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with the potential energy given by 



rrii^ij = —rriimj- — -. (2.26) 
For the interaction with cell k we write ()Goldstein 1980|1 

and 



1 1 m/c 
rrik^kj = -muruj^^^- -[Ykj -(^-Vkj) i;^^^^. (2.28) 

Here Q is the traceless quadrupole tensor about the centre of mass. It is defined as ([Goldstein 1980|1 

Qab= ^ rnp(3xa,pXb,p-rl6ab), (2.29) 
p=i 

where a, b run from 1 to 3 for each direction in Euclidean space and Sab is the Kronecker 
delta, e.g. Qu = Yl ^pC^x^ - Vp - z'^) and Q12 = Sj^^pXpyp, etc. 

The quadrupole moment of a cell in the tree is based on the quadrupole moments of its 
sub cells: 

^subcell -^subcell 

Qab= ^ iQab)p+ ^ mp{3Xa,pXb,p- rpab), (2.30) 
p=l p=l 

where Xa,p, Xb^p are now the coordinates of the subcells in the reference frame of the parent 
cell. 

With SPH we try to describe all quantities involved as continuous functions both in space 
and time. Therefore, we have to reduce the magnitude of close particle interactions in order to 
avoid unduly large gravitational accelerations (oc 1/r^). The particles are smoothed similarly 
to the way they are treated for hydrodynamics, i.e. as spherically symmetric and finite in 
extent with a radius of 2e. If two particles overlap, the mass involved in the calculation of the 
mutual gravitational acceleration is calculated from Gauss' gravitational theorem, otherwise 
the particles are treated as point masses. 

For particle i the mass interior to radius sei is given by m{s) = miW*{s) using pi{s) = 
mie~^W{s) and hence 
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W*{s) = [ 4Tru^W{u)du. 
Jo 



If particle i is less than 2eij = ej + ej away from particle j then the gravitational acceler- 
ation at particle j will be 

SigravAj = -miW* (2.31) 

which means that the mass of particle i outside r.y = |rjj| is not being taken into account. 
With this symmetric description of the gravitational interaction between particles i and j, 
we do not have to specify to which of the two particles we have applied Gauss' gravitational 
theorem. Then the total gravitational acceleration at particle j is 



E-'^Ml^ (2-32) 



The potential at distance rj, away from particle j is then given by 



r°° m ■ 

Jrij I 



Integrating by parts we obtain the expression for the mutual potential energy of two 
particles i and j: 



. = _!!!!!!^ iw*{s) + W**{s)) , (2.33) 

Tij 



where 



POO 

W**{s) = s 4TruW{u)du. 

J s 



For the M4-kernel we obtain 



WW) = ^ 



40s=^ - 36s^ + 15s^, O^s^l; 

80s^ - 90s^ + 36s^ -55^-2, 1 ^ s ^ 2; (2.34) 

30, s ^ 2, 
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and 



— < 

10 



s 



14 - 20s2 + I5s^ - 6s^ O^s^l; 
(2s + l)(2-s)^ 1 ^ s ^ 2; 



(2.35) 







s^2 



e can take either a small constant value (e- softening) which is not very good when the distance 
between particles changes considerably during the simulation; or it can take the value of the 
hydrodynamical smoothing length h, which has the benefit that it smoothes the gravitational 
and hydrodynamical forces by the same amount ({Bate &: Burkert 1997|) . We use the latter 
choice, i.e. = hij. 

2.6 Smoothing length 

In this section we shall discuss the importance of the hydrodynamical smoothing length, h, 
and we shall define the range of values that it should take. SPH is a Lagrangian particle 
method, with the particles - sample points moving with the fluid, containing constant mass 
within their radius of influence. In particular, the mass of the SPH particles is considered to 
be smoothed over a finite volume. As mentioned in ^2.21 at any time, only a few neighbouring 
particles overlap with the smoothing radius of any particle. Therefore, we choose the radius 
of influence of the smoothing kernel function, 2h, so that for any particle and at any time, 
it contains an approximately constant number of neighbouring particles. For this reason, we 
will use a time-varying non-universal h, as each particle requires its own smoothing length 
(e.g. in a dense region a smaller h is required, than in a rarefied region, in order to contain 
this constant number of neighbours). Nelson & Papaloizou p993t 11994(1 have shown that 
with adaptive h the energy is conserved quite accurately. 

Because we have to treat all the physical quantities involved as continuous functions 
both in time and space, we need to take into account a large number of interactions to 
reduce sampling errors. This needs to be balanced by the fact that for too big a smoothing 
length SPH will under-estimate the self-density of the particles and will over-smooth all the 
properties of the fiuid. Therefore the problem of finding the correct value for h is reduced 
to finding the right value for the number of neighbours, N^ant, for each particle. Tests on 
known distributions have shown that N^ant ~ 50 neighbours within 2h gives good results 
with the 3-D M4-kernel. 
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We allow a 10% fluctuation in the value of the number of neighbours, A^, for each particle 
at any time. Thus, we allow to be between Nmin = 45 and Nmax = 55. If N"^ is the number 
of neighbours of particle j at time-step n, then a trial value for hj at time-step (n + 1) is 
obtained by its value at the nth time-step according to 

hT' = h] (^^^ ' . (2.36) 

We then count the number of neig hbours within 2/i"+\ iV"+^ If Af"+^ is within the above 
limits then the value of ^j"*"^ is accepted. However, if the value of N'^^^ is not within the 
limits, we iterate over Eqn. 12.361 finding a new trial value for h using W-^^^ obtaining a 
new number of neighbours according to this new trial value of /i, and so on. We stop when 
is acceptable or when the fractional change of two successive trial values for h becomes 
less than 1%. After we have obtained the /i""*"^ values for all particles, we substitute these 
values to the current value of h. To set the initial values of h for all particles, we repeat the 
above procedure until either all particles have acceptable value for A^, or we exceed a finite 
number of iterations. Numerical tests indicate that after 20 iterations most particles have an 
acceptable A^. 

We will now discuss the way we identify neighbours. A direct search would be an A^^ 
procedure, but using the spatial information stored in the gravity tree we can speed it up 
considerably. The original idea ( Hernquist &: Katz 1989 1 for looking up in the tree for neigh- 



bours involved constructing a trial cube of side 4/ij for particle j and then identifying which 
cells this cube overlaps, always starting with the rootcell. If these cells contain subcells they 
are then examined in a hierarchical fashion until we open subcells containing single particles. 
We then compare the inter-particle separation with 2/ij and decide accordingly. 

We use a slightly modified method in order to avoid opening cells unnecessarily (jBhattal 1996|1 . 
We store in the tree the values of h for each particle contained in a cell. This way we can con- 
struct a bounding box, or 'kernel box', for each cell. This represents the minimum box that 
contains all the particles in the cell, when each particle is taken to extend to its smoothing 
radius hi (see Fig. 12. 2|) . The search process is as before but with a trial cube of 2/ij this time, 
using the kernel box for each cell instead of the cell itself. With this technique, we find all 
neighbours, missing no interactions within 2/ijj, Vi, j : i 7^ j, as hij is the value we are actually 
using for our SPH equations (Eqn. I2.19l fc r2.18|) . We only open a few cells unnecessarily and 
we can formulate this technique together with the calculation of the centres of mass for each 
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b) 



Figure 2.2. Example of the kernel box method for finding SPH neighbours. We use the minimum 
box that contains all the particles in the cell, when the particles are taken each to extend to its 
smoothing radius hi, instead of the cell itself (dashed box), a) Avoiding to open a cell unnecessarily: 
Particle j would interact with the cell but now it does not interact with the kernel box. b) Obtaining 
an interaction we would have missed: Particle j overlaps with particle i, but it does not overlap with 
the cell. 

cell so that we do not have to perform unnecessary walks of the tree. 

2.7 SPH equations 

We can now present the final forms of the three equations determining the evolution of the 
simulation. These are 



(2.37) 



dt 



E 



rrii- 



<-i] I 



+ 



2 + Hjj 



W 
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<-i] I 




(2.38) 



Pi 



1+1^) 

PoJ 



4/3- 



1/2 



(2.39) 



where Hij is given by Eqn. rTM For Eqn. rTM we have combined Eqns. EH ESI] & 12311 
W, W and W* are given by Eqns. 12.51 12.91 & 12.341 respectively. 
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Parenthetically, we note that if we were solving for u, the equation for the rate of change 
of the specific internal energy, after including heating from the artificial viscous forces, should 
read as 



-JT = n y^'^Ai \ ^ + ^ + ] '"ij-w x"^ r^+— y^^Ai ^ V- 

dt 2^ p] n \hijj\rij\ Pi \hi, 



2.8 Integration scheme 

We advance the positions and velocities of all particles in time using the second order Runge- 
Kutta integration scheme. This means that to advance particle j from the nth to the (n+l)th 
step, first we need to calculate its position and velocity at the midpoint. This is given by 



r 



n+l/2 



+ v] At/2 (2.40) 



^ (2.41) 

where is the total acceleration of particle j at the nth step and At is the discrete time-step 
with which all particles will be advanced (jPress et al. 199n|l . We calculate a" from Eqn. 12.381 
It is then easy to obtain the position of particle j at the (n + l)th step from 



rj^+i = r^" + ^2.42) 

However, in the meantime, we must calculate the total acceleration for j at the midpoint, 
since its velocity at the (n + l)th step is given by 

^n+l ^ + sT^+^l^ At. (2.43) 

The selection of the time-step At is of great importance. There are several time scales 
that can be defined locally in systems like the ones we follow. Firstly, the inverse of the local 
velocity divergence. Secondly, the ratio of the local length scale to the velocity at this scale. 
Thirdly, the square root of the ratio of the local length scale to the acceleration at this scale. 
And finally, a time scale similar to the second one, except that it involves the local sound 
speed instead of the total local velocity (separating the hydrodynamical properties of the 
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fluid from its overall behaviour). For each particle, i, we calculate the smallest of these time 
scales using its smoothing radius (radius of influence) as a local length scale, i.e. 



Ati = jMIN 



1 



hi_ 

IV • v\i ' \ vi 



(2.44) 



where 



(Ji = c, + C{aCi + p MAX {fi,j}) (2.45) 

j 

is a modified sound speed which includes the effect of artificial viscosity. C is a parameter 

usually taken equal to 1.2. a and /? are the viscosity parameters f ^2.4() . MAX gives 

j 

the largest contribution by a neighbour of particle i to its viscous acceleration. The value of 
7 is usually referred to as the Courant number and is given a sufficiently small value for the 
simulation to be well behaved as well as to conserve the total energy. We have adopted the 
value of 7 = 0.3. 

By choosing the smallest of these scales, we ensure that we do not evolve each particle for 
a time longer than any of the time scales dictated by the local dynamical/hydrodynamical 
properties. For the same reasons, we select the smallest of these minimum particle time scales 
for the value of the global time-step At at any step n. Formally this is given by 



TV 

At = MIN { 7 MIN 



i=l 



{ 



1 hi f hi_\^^^ 
IV • v\i ' \vi\ ' V |a,:| / ' a,- 



(2.46) 



2.9 Multiple time-steps 

If Eqn. I2.46l gives the value of the global time-step, it is obvious that there should be particles 
whose minimum local time scale is much larger than At. This means that if we can avoid 
evolving these particles with the minimum global time-step but with a time-step closer to the 
value they require, then there is a great gain in the speed of computation. This is the basic 
reason why we use the method of multiple time-steps (jBhattal 1996|) . This method is ideal for 
simulations where there exist within the computational domain both dense regions (requiring 
small time-steps) and rarefied regions (not requiring small time-steps). The particles are 
assigned individual time-steps which are allowed to vary from step to step according to their 
need. 
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Figure 2.3. Graphical representation of multiple time-steps for the example of TV = 5 time bins. The 
steps of particles in different time bins for a period of Atmax- A time-step At — Atmax/^ (time bin 
n = 2) is accepted for a particle, only at s = 0, 4, 8, 12. 

Taking into account the fact that at each step there are two half steps, the method 
creates a hierarchy of time-step bins, each containing particles that take one half step while 
the particles in the immediately lower bin take one full step. Therefore, the time-steps at any 
bin are twice as large as the ones at the immediately lower bin. The fact that the particles 
are not allowed to move with arbitrary time-steps but under this hierarchy of time bins, 
is dictated by the need to know the positions of all particles every time we calculate the 
accelerations and to therefore keep the system synchronised at regular intervals. 

The values of the discrete times-steps used by particles in different time bins, are therefore 
calculated as fractions of a maximum time-step, /Atmax- In particular, the time-steps can take 
the following values: Atmax ^ Atmax/'i, Atmax/^, Atmax/8, ... , Atmax f'^^''"'' ^- We choOSC 
the total number of available time bins, Nf^ins^ sufficiently large in order not to put any 
constraint on the time evolution of the simulation. Since at any time during a simulation 
/^tmax = Atmin'^^™'''"^^ , where Atmin IS the minimum time-step from the ones used at this 
time (corresponding to the A^mmth time bin), we can express the current position along the 
largest time-step as s Atmin, where s = 0, 1, 2, 3, . . . , 2^'""'^^. If s = the system is at the 
start of the maximum step and all particles are in-synch. 

For any particle, we calculate the ideal value of its time-step, AtideaU from Eqn. 12.441 
We then put this particle into the next smaller time bin, n, with time-step At = Atmax / '^'^ , 
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where n is defined as 

n = j^r| MAWAt.deaO |^^^ n = 0,l,2,3,...,(iV,,„,-l). (2.47) 

To ensure that at the end of Atmax all particles are in-synch, we have to check whether we 
can accept this bin, n, or not. A time bin is accepted only if the time from the beginning 
of the Atmax period is a multiple of the time this bin represents. Otherwise we choose the 
immediately lower acceptable time bin. For example, for N=5 we have Atmax = 16Atmm- 
If the time-step we are checking is At = Atmax/ ^ {n = 2), we can accept it only if s = 
0,4,8,12. Otherwise we will assign to this particle a time-step At = Atmax/8 {n = 3) if 
s = 0, 2, 4, 6, 8, 10, 12, 14, or the lowest available time-step (bin n = 4) if s is odd (Fig. 12. 3j) . 
With this test we ensure that a particle can move down to a lower time bin at any time, but 
can only move up the time ladder at times which allow the system to remain in-synch. 

As mentioned above, we need the positions and velocities of all particles to calculate the 
acceleration of the particles in the minimum time bin. To minimise errors we must update 
the positions and velocities of the particles in the upper bins for the intervening times. This 
is achieved by extrapolating the positions and velocities of these particles. Finally, at the 
end of every Atmax period the system is synchronised and the particles in all time bins have 
their positions and velocities updated by the integration scheme. 

2.10 Going through the code 'step by step' 

Fig. 12.41 shows a flow-chart of the algorithm that dictates the evolution of the fluid in time. 
It demonstrates the way we assemble all the previous features of the numerical method. It 
represents the cycle n of the integration scheme, when the system advances with the time-step 
in the minimum time bin, dt = Atmin- 

As mentioned in ^2.81 each step consists of two half steps. The first half step starts by 
advancing the system by Atmin/2. Active particles have their positions and velocities updated 
by the integration scheme (Eqns. 12.401 fc [2. 41|) . while all other particles have their positions 
and velocities updated by extrapolation. We then need to calculate the acceleration at the half 
time-step for the particles in the minimum time bin. We proceed as follows: we make a new 
tree since the particles have just moved ( H'2.b\i . Subsequently, we update the value of h for the 
particles in the minimum time bin f ^2.6() . We can then calculate the hydrodynamical f ^2.3() . 
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Figure 2.4. Flow-chart of the algorithm that dictates the evolution of the fluid in time. It represents 
the cycle, n, of the integration scheme, when the system advances with the time-step in the minimum 
time bin, dt = Atmin- 
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viscous (3231 and gravitational f ^2.5() acceleration for these particles. The total acceleration 
at the half time-step is given by Eqn. 12.381 

We then proceed to the second half step, when tasks similar to those of the first step 
are performed. Specifically, the second half step starts by advancing the system by another 
(Eqns. l2.42l fc E.43() . and thus completing a full time-step of the minimum bin Atmin- 
Once again, particles in the minimum time bin have their positions and velocities updated by 
the integration scheme (using the acceleration at the midpoint calculated during the first half 
step), while all other particles have their positions and velocities updated by extrapolation. 
We then find the active particles for the next time-step, i.e. the particles that are going to 
be updated by the integration scheme and therefore need their acceleration calculated. This 
time the active particles are not just those in the minimum time bin, but all the particles for 
which the total time of the simulation is a multiple of, either the whole time-step, or half the 
time-step of the time bin they are in. For example, for s = 4 in Fig. 12.31 all the particles 
in the three smaller time bins are active as they start a new step, as well as the particles 
in the second larger time bin as they have just completed their first half step. We will then 
calculate the acceleration for all the active particles, following exactly the same steps as in 
the first half step. Finally, before starting time-step (n + 1) we re-distribute the particles into 
time bins. In order to make sure that at the end of a Atmax period the system is in-synch, 
we allow particles to move up the time ladder only if Atmax is a multiple of the time-step of 
their time bin, which in our example means the particles in the three smaller time bins. All 
particles that have just completed a full step are allowed to move to smaller time bins, in our 
example again the particles in the three smaller time bins. 

Having calculated the accelerations for all the active particles, we can then proceed to 
the next time-step, (n + 1). 

2.11 Tests 

We shall now describe the performance of the above numerical method on some standard tests 
for a self-gravitating hydro dynamical code. In particular, first we simulate the interaction of 
two colliding flows in order to test the efficiency of the code's treatment to hydrodynamics. 
In particular, we would like to quantify the efficiency of artificial viscosity ( ^2.4|) . We then 
allow a uniform sphere to collapse freely to test the efficiency of the TCG method. Finally, 
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Figure 2.5. The density (left) and the velocity (right) of the simulated shock (points). The analytic 
solution is given by the solid lines (Dyson & Williams 1980). The shock is broader by less than 2h 
compared to the analytic solution. It is very well resolved (contains ~ 20 post-shock h). 

we combine both gravity and hydrodynamics to follow the evolution of a stable isothermal 
sphere in order to test the overall performance of the code. For each test we compare our 
results with the relevant analytic or semi-analytic solutions. We conclude that our numerical 
method treats effectively both the gravitational and hydrodynamical properties of the fluids 
we are going to simulate in chapters El & 13 

In quantifying the results of our tests, we need to be able to associate these results only 
with the performance of our numerical code. In order to decrease the numerical noise input 
by the initial distribution of particles, we perform all tests using clouds whose particles are 
taken initially to be on a lattice. 

2.11.1 Colliding Flows 

In order to demonstrate the efficiency of the code's treatment of hydrodynamics, we have 
chosen to simulate the interaction between two colliding flows instead of a Riemann shock 
tube. We believe that the former test is more relevant to the problems we will investigate 
later in this thesis, i.e. high Mach number isothermal shocks with high compression factors. 
The typical Riemann shock tube test only involves a compression factor of ^ 2. Nevertheless, 
Hosking p999|l has shown that the numerical method we are using, gives very good results 
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for the Riemann shock tube test. 

Two colHding isothermal flows will produce a strong shock, provided the Mach number of 
the collision, A4, is high. If uq is the pre-shock velocity of each flow in the reference frame of 
the shock and cq is the isothermal sound speed, then the Mach number is given hy A4 = uq/cq 
and we obtain the post-shock velocity of each flow, U2, from 

U2 = ^ (2.48) 

Uq 



(Dyson &: Williams 19801. From Eqn. 12.481 we obtain the value for the velocity of the shock. 



V^, with respect to the frame of the simulation, knowing the pre-shock and post-shock ve- 
locities of each flow in this frame of reference, vq and V2 respectively. Speciflcally, because 
= vq — Vs and U2 = —Vs (as V2 = 0), Eqn. 12.481 gives 



-V - 



vq - Vs 



s 



After some algebra we obtain 



V.= "±- M±±^, (2.49) 
From the continuity equation we can write 

92 = — 

n2 

where po 92 are the pre-shock and post-shock density for each flow, respectively. Finally, 
we can obtain an analytic estimate for the post-shock density for each flow 

92 = —r — Po- (2.50) 

Vs 

The simulation involves two colliding isothermal flows (T = lOK) each of unit length in 
the direction of the collision. Each flow has a velocity of 1 km s^^ in this direction. The 
isothermal sound speed at lOK is ~ 0.2 km s~^. Both flows have unit pre-shock densities. 
There are ~ 10, 000 particles in total, taken initially on a lattice. The other two dimensions 
have lengths of ~ Ah. Eqn. 12.491 gives Vg = 0.039 km s~^. Therefore, the Mach number of 
the shock is 7W ~ 5.20. This is typical of the values we shall later model. Eqn. 12.501 gives a 
post-shock density of p2 = 26.64. 

We evolve this collision with our 3-D SPH code, using periodic boundary conditions. We 
do not include the TCG method as it is not relevant for this test. Fig. 12.51 shows our results 
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at the point where the rarefaction waves at the opposite ends of the tube have propagated 
to one tenth of the initial length for each flow. The mean value for the post-shock density, 
P2, is ~ 26 having a ~ 11% dispersion around the mean. The mean value for the post-shock 
velocity, f2, is to the 3rd significant figure, with a dispersion of ~ 0.19 km s^^ around the 
mean. The shocked layer is not broader than the analytic solution. In fact, the wings of the 
shocked layer are extended to less than 2h. The layer is well resolved as its width is ~ 20 
times larger than the post-shock h. 

The simulation reproduces the analytic compression factor accurately. In fact, it is slightly 
lower than the analytic value, as some small fraction of the pre-shock kinetic energy is trans- 
formed to post-shock random velocity dispersion instead of work done in compression. How- 
ever, the biggest contribution to this random velocity dispersion comes from the breaking of 
the symmetry of the lattice. The particles in the shocked layer are not on a lattice any more 
and have random motions in all 3 dimensions. The artificial viscosity is very efficient though, 
as at the centre of the shocked layer the random motions are damped very effectively in all 
directions. Very little interpenetration is observed in the layer. 

We conclude that the results of this test are encouraging, and thus the choice of the values 
for the artificial viscosity parameters a = /3 = 1 is justified. The simulation has produced 
a well resolved, not broad shock with a compression factor very close to the analytic value. 
Comparison of this test with realistic simulations will follow in chapter 13 

2.11.2 Free-Fall collapse 

For the free-fall collapse simulation, we let a uniform density sphere of mass Mq and radius 
Rq collapse under its self-gravity. For this test we do not calculate the SPH equations, since 
the aim of this test is to verify the efficiency of the TCG method. Therefore, the acceleration 
in the integration scheme f ^2.8() is given only by Eqn. 12.321 

The sphere should collapse homologously to a point after a free-fall time, tff. The analytic 
solution for the free-fall collapse of a uniform sphere states that any fiuid element initially at 
radius ro will arrive at radius r at time t given by 





where ^ t ^ tff { Spitzer 1978 ). The free-fall time for a uniform density sphere is defined 
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Figure 2.6. Evolution of the 90%, 50% and 10% mass radii of a uniform sphere in free-faU collapse. 
The points show the values obtained from the simulation every 50 time-steps, while the solid lines 
indicate the analytic solutions. 

as 

We construct a Mq = 1 uniform density sphere with ~ 10, 000 SPH particles on 
a lattice. The sphere is then let to collapse for ~ 1 tjj. As particles come closer due to 
collapse, the gravity softening e = h is becoming smaller. However, h decreases slower for the 
particles close to the edge of the sphere than for those at the centre. This happens because 
particles close to the edge can find their ~ 50 neighbours only from the inner side of the outer 
boundary of the sphere. Having a larger h than they should, the gravitational acceleration for 
these particles is over-softened. This reduces the rate with which they collapse. To overcome 
this problem, and only for this test, we assign to all particles the same h = h, the mean value 
ofh. 

Fig. 12.61 compares the analytic solution (solid lines) with the values obtained from the 
simulations every 50 time-steps (points), for the 90%, 50% and 10% mass radii. The above 
edge effects, although reduced due to our treatment of h, give a small over-estimation for the 
90% mass radius. The inner two radii agree closely with the analytical solutions. 
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Figure 2.7. Initial radial density profile (left panel) and the radial density prohle after ~ 25 tff 
(right panel). The solution to Eqn. 12.541 is given by the solid lines. 

2.11.3 Stable Isothermal Sphere 

Finally, we test both the SPH and TCG methods in evolving a stable isothermal sphere. We 
construct the sphere using a combination of the hydrostatic balance and mass conservation 
equations 

1 d /r^dPX 

^^{-^)+i7rGp = (2.53) 
dr \ /O dr / 



with initial conditions 



piO)=Po, ^(0) = 0. 
dr 



Using the Chandrasekhar H1939|) dimensionless expression for Eqn. 12.531 we obtain 

«0)^|(0)^0, ,2.54) 

where we have substituted P = Cq p, p = poe~'^ and r = (47rG/9o)~^''^co^. The solution of 
Eqn. 12.541 formally extends to — > oo. We can truncate the sphere at a finite radius, H, 
provided that the pressure at the boundary of the sphere is balanced by an external pressure. 
For H < 6.45 the sphere is stable. 

We then solve Eqn. 12.541 numerically. We tabulate the values for ^, ip and dtp/dS,. Giving 
values to the total mass of the sphere Mq, the sound speed cq and H we can substitute back 
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to find r, p{r) and M(r). We have chosen Mq = 1 M©, cq = 0.17 km s ^ (corresponding to 
T = 7.9K) and E = 3. 

We then move ~ 10,000 particles taken from a uniform density lattice to reproduce the 
tabulated values for r, p(r) and M(r). The solution to Eqn. 12.541 is given by the solid lines 
in Fig. 12.71 We have also plotted the initial radial density profile (left panel) and the radial 
density profile after ~ 25 tjf (right panel) . We use crosses to present our initial data due to 
the fact that the initial distribution comes from a lattice and the radial profile consists of too 
few different points. We do not need crosses in the right panel of Fig. 12.71 as the particles 
have moved from their initial positions. However, the simulated radial density profile obeys 
the analytic solution very well. Note that at radii larger than ~ 0.04 pc, there is an under- 
estimation of the density due to boundary effects similar to those described in ^2.11.21 In 
particular, particles close to the edge of the cloud have larger h than they should. 
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Chapter 3 



Simulations of Rotating Clouds 
with m=2 Density Perturbations 

Simulations of fragmentation are only reliable if the Jeans condition is obeyed (Truelove et 
al. (flWTiriM) . Klein et al. (fl999|l . Boss et al. ^(2M^ . Bate & Burkert (fmoTll . Whitworth 
H1998|) : for a review see M.lf) . In this chapter we perform the standard test simulation first 
proposed by Boss & Bodenheimer (|1979|) . We show that our SPH code faithfully reproduces 
the results obtained by Truelove et al. 1)19971 119981) using an adaptive finite difference code, 
and by Bate & Burkert ()1997|1 using their SPH code ( ^3.2|) . This way, we test our code on a 
more realistic application and find that the results are consistent with those of Eulerian codes 
and other SPH codes. We also draw conclusions on the significance of the Jeans condition 
for fragmentation simulations. We note that the Bate & Burkert (|1997|1 SPH code has been 
developed independently from our code and differs from our code in several fundamental 
regards. 

We perform a series of simulations by changing the density above which adiabatic heating 
operates f ^3.4.2() . This way we can make a direct comparison between the results of our code 
and those of Bate & Burkert (fTMTjl . 

We also perform a series of simulations by gradually increasing the number of neighbours 
as a means of testing convergence to the known solutions ( ^3.4)1 . We perform both isothermal 
simulations and simulations which include adiabatic heating. For each case, there are 2 
subsets of simulations: low resolution and high resolution, depending on the total number of 
particles used. 

39 



40 CHAPTER 3. ROTATING CLOUDS WITH M=2 DENSITY PERTURBATIONS 



Knowing the solution to this problem, we can also explore some other numerical param- 
eters of the code. In particular, we verify the choice of the interpolating kernel (M4) by 
repeating a few simulations using a different choice of kernel ( ^3.5|) . We also test whether our 
results are corrupted by the choice of the initial spatial distribution of particles (i j.S.Bj) . But 
first, let us define the initial conditions for the simulations described in this chapter. 

3.1 Initial Conditions 

The standard test simulation for a Star Formation code is the evolution of a rotating, spher- 
ical, uniform-density, isothermal cloud with an m = 2 perturbation. The initial conditions 
we have used for this simulation are taken from Boss & Bodenheimer ()1979|) . In particular, 
we have constructed a uniform-density, isothermal spherical cloud of mass M = 1 Mq and 
radius R « 0.02 pc. The sound speed is cq = 0.17 km s~^ (corresponding to T = 7.9K), 
which assigns to the cloud a ratio of thermal to gravitational potential energy a ~ 0.26. The 
spherical cloud is constructed with particles cut either from a settled uniform distribution 
or from a uniform face-centred cubic lattice. The particles are then given an m = 2 az- 
imuthal perturbation with amplitude A = 10%, by adjusting their spherical polar azimuthal 
coordinate, (p, to a value (p* given by 

, , Asm{m(l)*) 

= H . 

m 

Finally, the cloud is given a uniform rotation (angular velocity O = 7.2 x 10^^^ rad s^^), 
so that the ratio of rotational energy to gravitational potential energy is /3 ~ 0.16. We use 
clouds having different numbers of particles depending on the problem. In all cases, we apply 
our full SPH code as given in chapter [21 When a cloud has to be evolved isothermally, we 
use 




instead of Eqn. 12.391 

In order to decrease the numerical noise input by the initial distribution of particles, we 
perform most simulations using clouds whose particles are taken initially to be on a lattice. 
To verify that the results of such simulations are not biased due to some preferred orientation 
of the initial lattice, we also perform one simulation where particles are taken initially from a 
"settled" distribution f p.6|) . Such a distribution of particles is produced when the particles 
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are taken in random positions and then they are relaxed to uniform density, using the SPH 
formulation described in Whitworth et al. ()1995|1 . 

All figures presented here are column density plots viewed along the rotation axis. The 
geometry of the problem (fragmentation happens on a flattened disc) supports the use of 
such plots, since projection effects are small. It also allows comparison with density contour 
plots or density equatorial slices, used by other workers, as most of the mass of the system 
ends up in the disc. The figure captions indicate the units of the colour tables. They also 
give the linear size of the figure and the time of the simulation. 

3.2 The solution in Eulerian and Lagrangian formulations 

The evolution of a rotating, spherical, uniform-density, isothermal cloud with an m = 2 
perturbation predicts that the cloud forms a flattened structure due to rotation, and that the 
inner region of the cloud initially expands; at all times, there are two over-dense zones due 
to the perturbation. After t ~1 tjf the inner region starts collapsing and forms an elongated 
high-density structure at t ~1.15 t//. At the two ends of the elongated structure material 
falls faster towards the over-dense regions. The two ends finally become self-gravitating at 
t ~1.20 tff and form a protostellar binary at t ~1.25 tff. The elongated structure between 
the binary components increases in density and forms a uniform density bar at t ~1.30 tff. 

The adaptive finite difference code of Truelove et al. (|1997t 11998(1 , implemented to obey 
the Jeans condition, shows that the bar does not fragment while the gas remains isothermal. 
In particular, Truelove et al. p998j) find that after t ~1.32 tjj the binary fragments are 
elongated and collapse to filamentary singularities (their Fig. 13), as suggested by Inutsuka 
&: Miyama (|1992() . Moreover, with their Fig. 12 they show that the bar between the binary 
does not fragment, contrary to what had been suggested by simulations using other grid codes 
(|Burkert ik. Bodenheimer 1993|) . They prove that fragmentation of the bar in these other 
simulations was a consequence of inadequate resolution. Therefore, since in their convergence 
tests, the resolution of their code can grow infinitely while resolving the local Jeans length 
and since with their highest resolution the bar does not fragment, they conclude that the bar 
should not fragment. Their results were subsequently confirmed by Boss et al. (I2UUU() . 

In Klein et al. (jl999j) the simulation is repeated and extended to higher densities with 
an equation of state that includes adiabatic heating with po = 5 x 10~^^ g cm^'^. Again, 
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the bar does not fragment, but due to adiabatic heating and thermal support the binary 
fragments now become spherical at t ~1.35 tff (their Fig. 2). They follow this binary for a 
few revolutions around the centre of the domain. The binary fragments accrete material from 
the bar and they finally form a detached binary at t ~1.45 tff (their Fig. 5). The remains 
of the bar form spiral arms around the rotating fragments (their Fig. 6). 

Bate & Burkert H1997|) repeat both the isothermal simulation and the simulation with adi- 
abatic heating using their SPH code. They, like Truelove et al. H1998|) . perform a convergence 
test. They increase the numerical resolution of their code by increasing the total number of 
particles in the simulations. In the isothermal regime, they obtain the expected evolution 
described above, with the simulation having the highest resolution (80,000 particles). They 
find that convergence appears to occur from the simulation with 40,000 particles. To prevent 
the simulation advancing with very small time-steps they use a minimum smoothing length of 
10^^ cm. Their isothermal simulations progress until t ~1.29 tff, where the Jeans condition 
is violated even for the simulation with the highest resolution (80,000 particles). Up to this 
point, the bar between the binary has not fragmented. They continue the simulation with the 
highest resolution (80,000 particles) including adiabatic heating^ initiating at po ~ 10~^^ g 
cm~^. This way the Jeans condition continues to be obeyed, as the Jeans mass increases with 
increasing density, due to the increase in temperature and thus in sound speed (c/. with Eqn. 
IA.8|) . They find that adiabatic heating provides the bar with extra support against collapse 
towards a filamentary singularity, and - in contrast to the finite difference simulations of 
Klein et al. (|1999j) - the bar fragments at t ~1.31 tff. This fragmentation is attributable to 
particle noise. 

They find that with higher pQ the bar becomes thinner and produces more fragments as 
the ratio of the length of the bar to its thickness increases, in accordance with the results 
of the grid code of Burkert & Bodenheimer H1993() . In particular, they show that fewer 
fragments are produced when more heating is applied, i.e. when adiabatic heating initiates 
at a lower density. They find that at t ~1.315 tff five fragments are produced when pQ = 
10~^^ g cm~'^, only one for = 3 x 10"^'* g cm~'^, and none for po = 10"^"^ g cm~'^. The 
survival or merger of the fragments depend on the chaotic dynamics of the system of these 
protostellar fragments. 

^They use a polytropic index of 5/2, implicitly including the two rotational degrees of freedom of H2. With 
their formulation, collapse is decelerated more slowly than with our code. 
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The figure is provided separately 
due to its large size 



Figure 3.1. Column density plot of the initial sphere using 80,000 particles. The linear size of this 
plot is 0.04 pc. The colour table has units of 1.18 x 10^ g cm^^. 

3.3 Changing po 

We have conducted the simulation for the evolution of a rotating, spherical, uniform-density, 
isothermal cloud with an m = 2 perturbation using our SPH code (321 • We have used a 
cloud of 80,000 particles with particles initially taken on lattice. We have included adiabatic 
heating. We have tried different values for po, the density at which heating starts, in an effort 
to explore whether convergence to the known solutions can be achieved with an equation of 
state that includes an adiabatic heating regime. Like Bate & Burkert ()1997() . we start with a 
value of Po = 10^^^ g cm~'^. We then repeat the simulations with smaller ( ^3.3.2)1 and larger 
f H3.3.3() values of po- We also set up a simulation using 600,000 particles that always obeys 
the Jeans condition and has the highest value for po = 5 x 10^^^ g cm~'^ ( ^3.3.4() . With this, 
we aim to use the least possible heating and let the simulation evolve isothermally as long as 
possible. Its high resolution enables us to draw conclusions on whether our previous results 
are resolution dependent. Table ITT] shows a summary of our findings. 

Fig. 13.11 shows the density projected on the x-y plane initially (at t = 0), where the 
density enhancements indicate the m=2 perturbation. 
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The figure is provided separately 
due to its large size 



Figure 3.2. Column density plots for a cloud of 80,000 particles before heating starts and at the end 
(po = 10~^^ g cm~^). The linear size of these plots is 0.004 pc. The colour table has units of 1.18 x 
10^ g cm~^. Top: Column density plot before heating starts {t = 1.25 tff). Bottom : Column density 
plot at the end {t = 1.30 tff). 
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The figure is provided separately 
due to its large size 



Figure 3.3. Column density plot for a cloud of 80,000 particles at the end {po = 5 x 10^^^ g cni^''). 
The time is t = 1.31 tff. The linear size of this plot is 0.004 pc. The colour table has units of 1.18 x 
10^ g cm~2. 



The figure is provided separately 
due to its large size 



Figure 3.4. Column density plot for a cloud of 80,000 particles at the end {po = 10~^^ g cm~^). The 
time is t = 1.33 tff. The linear size of this plot is 0.004 pc. The colour table has units of 1.18 x 10^ 
g cm~^. 
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3.3.1 po = 10"^^ g cm^3 

The results are similar to those of Bate & Burkert (|1997|1 . The simulation evolves isothermally 
until t ~1.25 tff producing a binary and a thin bar. At this stage, the binary components 
are elongated (top panel of Fig. \'A.2\i . 

After this point adiabatic heating initiates. Subsequent collapse onto the binary com- 
ponents is decelerated. Each of these objects has thermal support. They are forced into 
solid-body rotation by the high effective shear viscosity. As the rotating elongated objects 
interact with the slowly infalling envelope, they form larger prolate objects with thermal 
support. At a later stage, the binary components become spherical. The spherical objects 
are followed by thin spiral tails. At t ~1.28 tff, the bar - which also has thermal support — 
fragments. At the end of the simulation (bottom panel of Fig. 13.2(1 . at t ~1.30 tff, the bar 
has produced 3 fragments, one at the bar's centre and two closer to the binary components. 
The binary components are likely to merge with the bar fragments. There are also 2 smaller 
fragments in the spiral tails, one in each. The mass for each of the binary components is 0.1 
Mq and the radius is 39 AU. Their separation is 330 AU. The total mass of the bar fragments 
is 0.02 Mq. The simulation has reached a peak density of Ppeak = 2.1 x 10^^^ g cm^'^. 

The simulation could resolve fragmentation up to a density of pmax = 9.6 x 10~^^ g 
cm^^. Therefore, fragmentation may be slightly unresolved in this simulation (see discussion 
in ^4.1() . We have repeated the above simulation with adiabatic heating initiating at po = 
9 X lO"^'' g cm~^. We obtain exactly the same results as above, with the bar fragmenting 
again at t ~1.28 tff. Therefore, we conclude that the bar fragments produced above are not 
due to the Jeans condition not being obeyed. 

We have repeated the simulation after changing pQ, the density where adiabatic heating 
initiates. We have used both smaller [po = 5 x 10~^^ g cm~'^ and po = 10~^^ g cm~^) and 
larger (po = 5 x 10"^'^ g cm~^) values than above. 

3.3.2 Decreasing po 

The simulation with po = 5 x 10~^^ g cm~^ produces the expected results (c/. the corre- 
sponding simulation of Bate & Burkert (|1997)1 with po = 3 x 10~^^ g cm~^). Since more 
heating is applied than before, the bar fragments later, at t ~1.285 tff, as it becomes thin 
later. At the end of the simulation (Fig. 13. 3|) . at t ~1.31 tff, the peak density is lower than 
before, Ppeak = 9-5 x 10^^^ g cm^'^, due to the increased amount of heating. The bar has 
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fragmented to only one fragment at the centre. The simulation has advanced further in time 
and the binary components have approached closer. There are some lumps in the spiral tails. 
The binary components each have mass 0.13 Mq and radius 33 AU. Their separation is 270 
AU. The mass of the bar fragment is 0.03 Mq. 

In the simulation with po = 10~^^ g cm~^ the bar does not fragment. The binary fragments 
have thermal support from an earlier time and thus become thicker. At the end of the 
simulation, t ~1.33 tff (Fig. 13. 4|) . the peak density is even lower than before, Ppeak = 1-8 
X 10~^^ g cm^^. The binary fragments have merged and completed 3/4 of a full rotation. 
The merging of the fragments is probably due to the artificial shear viscosity being too large 
with our implementation of viscosity ( ^2.41 - see discussion in ^3.3.5|1 . This way, we cannot 
reproduce the binary system followed by Klein et al. p999j) . 

3.3.3 Increasing po 

For the simulation with the least heating (po = 5 x 10~^^ g cm~^) the results change in the 
opposite sense: the bar fragments earlier {t ~1.270 tjf) and at the end of the simulation, t 
~1.275 tff (Fig. 13. the bar has broken into 9 fragments; there are also 2 lumps in each 
spiral tail. The masses of the binary components are 0.05 M0 and their radii are 19 AU. 
Their separation is 515 AU. The total mass of all bar fragments is 0.06 Mq. We cannot rule 
out artificial fragmentation for this simulation as we did above, since at po = 5 x 10~^^ g 
cm~'^, the Jeans mass has stopped being resolved for half a decade in density. At t ~ 1.275 
tff, the peak density has reached Ppeak = ^-7 x 10^^'^ g cm~^, a value much higher than 
before. 

To resolve fragmentation up to this high density (po = 5 x 10'^^ g cm~^) we need 185,000 
particles (c/. Eqn. 14.3)1 . We have conducted such a simulation and the results are similar 
to those of the simulation with po = 10"^'^ g cm~^. In particular, the binary and a thin 
bar formed. The mass of each binary component is 0.04 Mq and its radius is 12 AU. Their 
separation is 495 AU. At t ~1.267 tff the bar fragments. There are 7 fragments in the bar 
plus one lump in each spiral tail, at the end of the simulation (t ~1.271 tff - Fig. 13. 6() . 
The total mass of all bar fragments is 0.02 Mq. Since the Jeans mass is always resolved, 
fragmentation in this simulation is realistically modelled {ppeak = 1-5 x lO"^*' g cm~'^). 
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The figure is provided separately 
due to its large size 



Figure 3.5. Column density plot for a cloud of 80,000 particles at the end (po = 5 x 10^^"^ g cm^'^). 
The time is t — 1.275 t//. The linear size of this plot is 0.004 pc. The colour table has units of 1.18 
X 10^ g cm~^. 

3.3.4 600,000 particle simulation 

Finally, we have conducted a simulation with 600,000 particles where adiabatic heating initi- 
ates at /?o = 5 X lO"^'^ g cm~^ (we use such a high density for the switch to adiabatic heating 
in order to evolve the cloud isothermally as long as possible). At the end of the simulation 
(t ~1.245 tff) the bar has not fragmented. The simulation has reached a peak density of 
Ppeak = 1.8 X 10~^ g cm~^, the highest of all the simulations we conducted. As Truelove 
et al. H1997|l suggest, a coarse simulation is less evolved at the same time compared to a 
fine simulation. This happens in SPH as well, since with higher resolution the particle h is 
smaller and thus the density modelled becomes higher. Truelove et al even suggest that 
between simulations of different resolution comparison should be made when they both have 
advanced to the same density and not to the same time. 

The bottom panel of Fig. 13.71 is a column density plot at the end of the simulation (t 
~1.245 tff). The bar is very sparse and it has just started becoming self-gravitating. Until a 
few time-steps before, the binary components were still rather elongated apart from their very 
centres where a spherical core had developed. This core is spherical mainly due to the fact 
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The figure is provided separately 
due to its large size 



Figure 3.6. Column density plot for a cloud of 185,000 particles at the end {po = 5 x 10~^^ g cni~^). 
The time is t — 1.271 tff. The linear size of this plot is 0.004 pc. The colour table has units of 1.18 
X 10® g cm~^. 

that the thickness of this elongated self-gravitating object is smaller than the /i of a particle 
in it (these are the particles with the smallest h in the whole simulation being the particles 
with the highest density). This way, due to the spherical symmetry of the kernel, the centre 
of each elongated object settles to form a spherical object containing its ~50 neighbours. 
The high effective shear viscosity has since then put the binary components into solid-body 
rotation and they have grown in size due to interaction with the accretion flows. At the end, 
the masses of the binary components are 0.008 Mq and their radii are 3 AU. Their separation 
is 680 AU. 

The bottom panel of Fig. 13.71 compares well with Fig. 2 of Klein et al. p999)l . In fact, 
our simulation has reached a higher peak density as it has evolved isothermally for 2 orders 
of magnitude more than that of Klein et al. The top panel of Fig. 13.71 is a column density 
plot just before adiabatic heating starts {t ^-1.237 tff). It compares very well with Fig. 12 
of Truelove et al. H1998|) . Again the peak density {ppeak = 5-2 x 10^^^ g cm~'^) has reached 
a higher value than that of the simulation of Truelove et al., as adiabatic heating in our 
simulation starts one order of magnitude higher in density than in theirs. 
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Figure 3.7. Column density plots for a cloud of 600,000 particles before heating starts and at the 
end (po = 5 X 10~^^ g cm~^). The linear size of this plot is 0.004 pc. The colour table has units of 
1.18 X 10^ g cm~^. Top: Column density plot before heating starts {t = 1.237 tff). Bottom : Column 
density plot at the end {t = 1.245 tff). 
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Table 3.1. Summary of results for the simulations with different po- For each simulation the third 
column gives the time of the binary formation, the fourth the final time of the simulation and the 
sixth the time of the fragmentation of the bar. All times are quoted in units of the free-fall time of the 
initial cloud, tff. The fifth column gives the number of bar fragments at the end of each simulation. 
The last column gives the peak density at the end of each simulation. 

3.3.5 Conclusions 

In general, our SPH code faithfully reproduces the results of the 80,000 particle simulations 
of Bate & Burkert (|1997)1 . We, hke they, conclude that with higher the bar becomes 
thinner and produces more fragments as the ratio of the length of the bar to its thickness 
increases, i.e. we show that fewer fragments are produced when more heating is applied. In 
addition, with the 600,000 particle simulation our code faithfully reproduces the results of 
Truelove et al. (|19971 11998|l and Klein et al. p999|) . This makes us confident that our code 
models self-gravitating gas dynamics realistically and therefore, we can use it in our effort to 
implement a method for obtaining higher resolution (chapter^, as well as in simulations of 
clump-clump collisions (chapter [nj. Table ITT] shows a summary of results for the simulations 
with different po- 

From this series of simulations, there are several conclusions we can draw. Firstly, the 
fragmentation of the bar in the simulations of Bate & Burkert (|1997|1 and our coarse sim- 
ulations (80,000 particles) is due to poor sampling. The fact that Bate & Burkert p997|l 
find that with increasing resolution the bar produces more fragments, cannot prove that the 
bar will always fragment, as they have not achieved convergence with their 80,000 particle 
simulation. Our 600,000 particle simulation suggests the opposite, i.e. the bar should not 
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fragment, in accordance with Truelove et al. p9971 [TI)98() . Unfortunately, we cannot predict 
bar fragmentation theoretically in a manner similar to the prediction of the fragmentation 
of a layer ()Whitworth et al. 1994al IWhitworth et al. 1994b|l . as it is not possible to make a 
dimensional analysis in one dimension. Truelove et al. (fTWzl [T998jl and Klein et al. (TTHMjl 
with their adaptive finite difference code can regulate their grid size achieving theoretically 
infinite resolution, a feature that SPH codes do not possess. One of the primary aims of the 
PhD project discussed in this thesis is to invent a method that would allow the resolution of 
an SPH simulation to be increased on-the-fly by particle splitting (see chapter 0J). 

Secondly, there may be a way of preventing bar fragmentation even for simulations with 
80,000 particles. This could be done with the implementation of a triaxial kernel, where h is 
replaced by a tensor that gives different smoothing lengths in different directions. The values 
for h in different directions are such that a particle still overlaps with ~50 of its neighbours 
(|()wen et al. 19981) . 

Thirdly, the times for bar fragmentation are based on our perception of a bar fragment 
to be a density enhancement on the bar twice as large as the underlying bar density. It has 
been shown that in all cases, such an enhancement will form a fragment, due to the elongated 
geometry of the bar and the spherical symmetry imposed by the kernel. This is an empirical 
law for this kind of simulation. We have considered the time that such a density enhancement 
needs to grow as the time of bar fragmentation. Our definition may differ from the relevant 
definition of Bate & Burkert (|1997j) and this may be the reason for the small mismatch in 
the quoted bar fragmentation times. 

Fourthly, the fact that the fragments of the bar merge with the binary components in all of 
the simulations that the bar fragmented, is probably due to excess shear viscosity introduced 
by our application of artificial viscosity f ^2.4|) . This is also the reason for the merger of the 
binary components in the simulation with po = 10~^^ g cm^'^. In particular, our code has 
failed the test for the collapse simulation of a rotating unperturbed sphere, suggested by 
Norman, Wilson & Barton p980|) . Shear viscosity acts as an agent that redistributes angular 
momentum and we could not obtain a perfect singular disc. Instead a ring was produced 
around the disc centre. 

An alternative implementation could be explored. One could use the switch introduced by 
Morris & Monaghan (|1997j) . We refer to their Eqn. 30 with the Balsara (|1995j) source term 
- Eqn. 4. Another switch that calculates viscosity only for physically approaching streams 
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could also be used, e.g. using V(V • v). 

Finally, it is interesting to note that bar fragmentation happens in a symmetric fashion, 
with a central object and equal number of fragments on either side of it, at equal distances. 
This can be attributed to the symmetric initial distribution of particles (taken from a lattice) 
that has removed some of the numerical noise and prevented the randomness in the fragment 
positions present in the simulations of Bate & Burkert (|1997|) and ^3.61 

3.4 Changing the number of SPH neighbours 

In ^2.61 and all simulations so far, we have used a fixed value for the number of neighbours, 
Nn ~50, contained within the smoothing radius, h, of all particles. The choice of this value 
for Nn was dictated by the need to balance accuracy - that requires large Nn ~ against 
computational cost which is reduced with small values for Nn- We have shown that medium 
resolution (80,000 particles) SPH implementations of the standard test simulation cannot 
produce the expected evolution for the bar between the binary components and that we need 
to increase the number of particles, N, by one order of magnitude before the results converge 
to those of the finite difference code of Truelove et al. (11997(1 . We have tried increasing not 
just the number of particles, N, but also the number of neighbours, N^, this is suggested by 
Rasio (|1999|) as a means of increasing SPH accuracy (reducing sampling noise). We use the 
limiting case of Rasio's suggestion by increasing A^^^ and with the same rate, in an effort to 
keep the same resolution (constant N/Nn) for all simulations and therefore make comparison 
between them more meaningful. 

We have conducted a series of simulations by increasing A'^ from 50 to 200 in steps of 50. 
N has changed respectively from 80,000 to 320,000 in steps of 80,000. In order to be able 
to identify the exact point where our simulations converge to the expected solution, we have 
also set up a series of low resolution simulations with N ranging from 10,000 to 40,000 in 
steps of 10,000. All these simulations are evolved both with an isothermal and an adiabatic 
equation of state. 

We present each of the four groups of simulations separately in order to avoid confusion 
and to identify clearly the effect of increasing A'^^ on the final state of the simulations. For 
consistency, we present column density plots for all simulations, viewed along the rotation 
axis. All plots have the same linear size as the plots of ^3.31 
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N 


Nn 


hin 1 tff 


tend 1 tff 


Bar fragments 


tfrag / tff 


Ppeak 1 g Cm-3 


10,000 


50 ±5 


1.265 


1.274 


2 


1.272 


3.2 X 10-"^ 


20,000 


100 ±10 


1.260 


1.278 





N/A 


2.5 X 10-8 


30,000 


150 ±15 


1.259 


1.276 





N/A 


6.0 X 10-9 


40,000 


200 ±20 


1.258 


1.274 





N/A 


2.7 X 10-9 



Table 3.2. Summary of results for the low-resolution isothermal simulations with increasing N^- For 
each simulation the third column gives the time of the binary formation, the fourth the final time of 
the simulation and the sixth the time of the fragmentation of the bar. All times are quoted in units of 
the free-fall time of the initial cloud, tff. The last column gives the peak density at the end of each 
simulation. 

3.4.1 Isothermal Simulations 

3.4.1.1 Low resolution isothermal simulations 

A summary of the low resolution isothermal simulations is given in Table l3^ All simulations 
can resolve fragmentation up to a density of pmax = 1-75 x 10-^^ g cm"'^. They are all 
terminated when the multiple time-step method ( H2.9|) runs out of time bins. In fact, all 
simulations stop when the time-step becomes less than 2 x 10-^ t//, so that it would have 
been computationally inefficient to continue (i.e. in order for time to progress by 10"^ tjj 
we would need 5-6 times the run-time up to that point). The four simulations have shown 
the following: 

1. A^=10,000 and A'^„=50 Due to poor sampling, both the initial expansion phase and 
the collapse that follows it are not properly modelled. This causes the binary to form 
later than expected, at t\,in ~1.265 tff (cf. the values of the third column of Table m|) . 
In fact, a filament forms first and both ends of the filament become self-gravitating 
shortly after, forming the binary (in accordance with the evolution of the 10,000 particle 
simulation of Bate & Burkert (|1997|1 ). The bar fragments shortly after the binary 
formation, at tfrag ~1.272 tff. The top left hand panel of Fig. 13.81 is a column density 
plot at the end of the simulation, tgnd ~l-274 tff. The binary components appear to be 
elongated. The peak density of the simulation has reached a non-physical value, Ppeak 
= 3.2 X 10^'' g cm^^. This simulation does not reproduce the expected evolution. 
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Figure 3.8. Column density plots for the low-resolution isothermal simulations with increasing Nn- 
Final stage of the simulations with 7V=10,000 and A'^;^=50 (left), iV=20,000 and iV„ = 100 (right) on 
the top row and iV=30,000 and Nn^lbO (left), iV=40,000 and A'^;^=200 (right) on the bottom row. 
The details for each simulation are given in Table 13.21 The results converge after the simulation with 
A^=30,000 and A^„=150. The linear size of all plots is 0.004 pc. The colour table has units of 1.18 x 
10^ g cm~^. 
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Figure 3.9. Column density plots for the high-resolution isothermal simulations with increasing Nn- 
Final stage of the simulations with iV=80,000 and N„=50 (left), A^=160,000 and 7Y„=100 (right) on 
the top row and 7V=240,000 and 7V„=150 (left), 7V=320,000 and A'^;^=200 (right) on the bottom row. 
The details for each simulation are given in Table 1531 The results converge after the simulation with 
A^=160,000 and A':„=100. The hnear size of aU plots is 0.004 pc. The colour table has units of 1.18 x 
10^ g cm"2. 
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2. A^=20,000 and Nn=100 The overall evolution of the cloud is similar to the that of the 
previous simulation. However, the binary forms a bit earlier than before and the binary 
components start collapsing almost simultaneously with the bar and not after the bar. 
At the end of the simulation, tend ~1.278 tff (top right hand panel of Fig. 13. 8|) . the bar 
has not fragmented although it is very dense and it might have fragmented if we had 
continued the simulation. At this stage, the binary components are rather elongated 
apart from a spherical core that has developed at their centres. The reason for the 
formation of such a core is discussed in ^3.3.41 The peak density at the end, Ppeak = 
2.5 X 10~^ g cm~'^ is lower than before, therefore the density field is more realistically 
modelled, as expected for a simulation with less noise (more particles contribute to each 
SPH sum and all quantities are better modelled). 

3. A^=30,000 and A^„=150 The binary forms even earlier and its components end up 
elongated with spherical cores formed at their centres. The bar is not so dense as 
before and clearly it has not fragmented by the end of the simulation, at tend ~1.276 
tff (bottom left hand panel of Fig. 13.8(1 . The peak density, Ppeak = 6.0 x 10~^ g cm~^, 
is even lower. 

4. A'^=40,000 and iV„=200 The binary still forms a bit earlier and its components have 
the same shape as in the previous simulation. The peak density at the end, Ppeak = 
2.7 X 10^^ g cm^^, is even lower. The bar does not fragment and it is even less dense. 
The bottom right hand panel of Fig. 13.81 is a column density plot at the end of the 
simulation, at tend ~l-274 tff. The evolution of this simulation is not much different 
from that of the previous simulation. Therefore, we believe that convergence is achieved 
from the simulation with A^=30,000 and A'„=150. 

Obviously, this set of low resolution simulations cannot reproduce the exact solution. It 
appears that low resolution in the early stages delays binary formation. 

We have also conducted a simulation with A^=40,000 and Nn=50 (having resolution four 
times higher than that of the above simulations) in order to confirm that the above results 
depend not just on the increasing number of particles but also on the increasing number of 
neighbours. With A^=40,000 and A',j=50, we have found that the binary forms earlier, at 
tbin ~ 1.252 tff, in accordance with the corresponding simulation of Bate & Burkert (jl997j) . 
The peak density at the end of the simulation with A'"=40,000 and Nn=50 is two orders of 
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N 


Nn 


hin 1 tff 


tend 1 tff 


Ppeak 1 g cm-3 


80,000 


50 ±5 


1.248 


1.252 


1.7 X 10"*^ 


160,000 


100 ±10 


1.247 


1.250 


3.7 X 10-9 


240,000 


150 ±15 


1.247 


1.250 


7.3 X 10-*^ 


320,000 


200 ±20 


1.247 


1.250 


8.8 X 10"*^ 



Table 3.3. Summary of results for the high-resolutfon isothermal simulations with increasing N^. 
For each simulation the third column gives the time of the binary formation and the fourth the final 
time of the simulation. All times are quoted in units of the free-fall time of the initial cloud, t//. The 
last column gives the peak density at the end of each simulation. 

magnitude higher than that of the A^=40,000 and A'^„=200 simulation. This is partly due to 
the higher resolution (a finer simulation is more evolved) and partly due to the increased noise 
(less particles contribute to the SPH sums) of the A^=40,000 and A'^„=50 simulation. The 
difference in tun and Ppeak in the A^=40,000 and A''„=50 simulation confirms the dependence 
of the results of the above four simulations on the increasing number of neighbours. 

3.4.1.2 High resolution isothermal simulations 

A summary of the high resolution isothermal simulations is given in Table HOI All simulations 
can resolve fragmentation up to a density of Pmax = 9.6 x 10^^^ g cm~'^. They are all 
terminated when the multiple time-step method ( ^2.91) runs out of time bins. In fact, all 
simulations stop when the time-step becomes less than 10^^ iff, so that it would have been 
computationally inefficient to continue (i.e. in order for time to progress by 10^^ iff we 
would need 8-10 times the run-time up to that point). 

All four simulations have shown that the binary forms at the expected time (tbm ~1.25 
tff) and that its components are elongated apart from their centres where a spherical core 
develops due to the spherical symmetry of the kernel (see discussion in ^3.3.4(1 . The binary 
components become thinner with increasing A^^, so that the simulations converge towards 
the expected solution of filamentary singularities with the noise being reduced. 

The bar between the binary components does not fragment, is sparse and becomes less 
dense with increasing Nn- The peak density at the end of the simulations becomes lower with 
increasing A^^^, as the density field is more smoothed with less noise. Note that the last two 
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simulations initially ended a bit earlier than the quoted values of tend in Table 13.31 We have 
restarted them by resetting the time bin hierarchy and continued them for a few time-steps. 
The binary components are in free-fall collapse so that the peak density has increased by 
more than one order of magnitude during this short time. 

Fig. 13.91 shows the end state of all four simulations. Convergence to the expected so- 
lution starts to appear from the simulation with A^=160,000 and A'„=100. This is a clear 
computational gain as we only need to double and Nn with respect to the 80,000 particle 
simulations of ^3.31 to obtain a solution similar to that of the finite difference simulations of 
Truelove et al. (|1997t I1998|l . Comparison with the low resolution isothermal simulations has 
shown that the increased resolution (8 times higher) has a clear effect on the morphology of 
all structures formed; the binary components appear more well-defined and the bar is more 
sparse. In the series of high resolution simulations convergence is achieved with a lower value 
for Nn. 

However, following the simulations at such high densities with an isothermal equation of 
state is not physical and in order to model the evolution of the cloud realistically we need to 
include adiabatic heating. This is presented in the following subsection. 

3.4.2 Simulations with Adiabatic Heating 

3.4.2.1 Low resolution simulations with adiabatic heating 

A summary of the low resolution simulations with adiabatic heating is given in Table 13.41 
All simulations can resolve fragmentation up to a density of pmax = 1-75 x 10~^^ g cm~^. 
Adiabatic heating starts at /jq = 5 x 10^^"^ g cm^'^. All simulations are terminated at tend 
~1.29 tff (Fig. I3l0l) . 

All four simulations have shown that the binary forms later than expected (tbm ~1.26 
tff) due to the low resolution that leads to inadequate modelling of the initial stages of the 
cloud evolution. After heating switches on, the binary components obtain thermal support 
and become spherical. The excess shear viscosity puts them in solid-body rotation and 
they eventually grow in size and come closer together (see discussion in H3.3.5I) . The binary 
parameters are similar to those of the simulations of ^3.31 

The whole cloud obtains thermal support and the peak densities do not reach such high 
values as in the isothermal simulations. The peak density decreases with increasing Nn. The 
last two simulations have progressed a bit more in time, therefore their peak density at the 
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N 


Nn 


hin 1 tff 


tend 1 tff 


Bar fragments 


tfrag / tff 


Ppeak 1 g 


10,000 


50 ±5 


1.265 


1.295 


4 


1.279 


5.7 X 10-1° 


20,000 


100 ±10 


1.260 


1.293 


9 


1.281 


4.2 X 10-1° 


30,000 


150 ±15 


1.259 


1.297 


9 


1.283 


4.7 X 10-1° 


40,000 


200 ±20 


1.258 


1.294 


9 


1.284 


4.6 X 10-1° 



Table 3.4. Summary of results for the low-resolution simulations with adiabatic heating (po = 5 x 
IQ-^'^ g cm^'^) and increasing A^„. For each simulation the third cohunn gives the time of the binary 
formation, the fourth the final time of the simulation and the sixth the time of the fragmentation of 
the bar. All times are quoted in units of the free-fall time of the initial cloud, tff. The last column 
gives the peak density at the end of each simulation. 

final stage is a bit higher than that of the simulation with A^=20,000 and A'„=100. 

In all simulations the bar fragments. It fragments at a later stage with less noise. However, 
decreasing the noise does not decrease the number of bar fragments as one might expect. 
The simulation with A^=10,000 and A'„=50 produces 4 bar fragments but the other three 
simulations produce 9 fragments. Some of these fragments merge with the binary components. 
In all four simulations a few lumps are produced in both spiral tails. 

Bate & Burkert H1997() have argued that with increasing the bar produces more frag- 
ments as it becomes thinner and the ratio of its length to its thickness increases. One could 
argue the same for increasing A^^. Moreover, since the Jeans condition is not obeyed in all 
four simulations, lumps containing less than N^ particles can start condensing out and it is 
easier to find such lumps with increasing Nn- 

We conclude that from the simulation with A^=20,000 and A^„=100 convergence is achieved, 
but clearly not to the right solution as given in Klein et al. (1999). 

3.4.2.2 High resolution simulations with adiabatic heating 

A summary of the high resolution simulations with adiabatic heating is given in Table ESI 
All simulations can resolve fragmentation up to a density of pmax = 9.6 x lO-i^ g cm-^. 
Adiabatic heating starts at po = 5 x 10"!^ g cm-^ in order for the Jeans condition to be 
obeyed. All simulations are terminated at tend ~1.31 tff (both top row panels and bottom 
left hand panel of Fig. I3.11j) . 
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Figure 3.10. Column density plots for the low-resolution simulations with adiabatic heating (po = 
5 X 10^^^ g cm"'^) and increasing iV„. Final stage of the simulations with A^=10,000 and A^„=50 
(left), 7V=20,000 and 7V„=100 (right) on the top row and iV=30,000 and iV„=150 (left), iV==40,000 
and iV„=200 (right) on the bottom row. The details for each simulation are given in Table The 
results converge after the simulation with Af=20,000 and iV„=100. The linear size of all plots is 0.004 
pc. The colour table has units of 1.18 x 10^ g cm~^. 
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Figure 3.11. Column density plots for the high-resolution simulations with adiabatic heating {po = 5 
X 10^^''' g cm^'^) and increasing Nn. Final stage of the simulations with Af=80,000 and Af„=50 (left), 
7V=160,000 and 7V„=100 (right) on the top row and iV=240,000 and 7V„=150 left on the bottom row. 
The details for each simulation are given in Table The results converge after the simulation with 
-/V=160,000 and A^„=100. The right panel of the bottom row is a column density plot for the high 
resolution simulation with Af=160,000 and 7V„=100 and delayed adiabatic heating (po = 5 x 10^^^ g 
cm"-^). The linear size of all plots is 0.004 pc. The colour table has units of 1.18 x 10^ g cm~^. 
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N 


Nn 


hin 1 tff 


tend 1 tff 


Bar fragments 


tfrag / tff 


Ppeak 1 g Cm-3 


80,000 


50 ±5 


1.248 


1.307 


1 


1.285 


9.5 X 10-12 


160,000 


100 ±10 


1.247 


1.316 


1 


1.294 


1.2 X 10-11 


240,000 


150 ±15 


1.247 


1.309 


1 


1.300 


1.2 X 10-11 



Table 3.5. Summary of results for the high-resolution simulations with adiabatic heating (po = 5 x 
IQ-i^ g cm-'^) and increasing Nn- For each simulation the third column gives the time of the binary 
formation, the fourth the final time of the simulation and the sixth the time of the fragmentation of 
the bar. All times are quoted in units of the free-fall time of the initial cloud, i//. The last column 
gives the peak density at the end of each simulation. 



The binary forms at the expected time {tun ~1.25 t//) and its physical parameters are 
similar to those of the simulations of ^3.31 After heating switches on, the binary components 
obtain thermal support and become spherical. The excess shear viscosity puts them in solid- 
body rotation and they eventually grow in size and come closer together (see discussion in 

The whole cloud obtains thermal support and the peak densities do not reach such high 
values as in the isothermal simulations. The peak density decreases with increasing N^. The 
last two simulations have progressed a bit more in time, therefore, their peak density at the 
final stage is a bit higher than that of the first simulation. 

In all simulations the bar produces one fragment at the centre. The bar fragments at a 
later stage with less noise. In the simulation with A^=80,000 and A'„=50 each spiral arm 
produces a small lump. There are no lumps in the spiral tails in the other two simulations. 
Therefore, we conclude that convergence is achieved from the simulation with A^=160,000 and 
A^n=100. Since convergence is achieved, we have not conducted a simulation with A^=320,000 
and A^„=200 to avoid the most computationally expensive of this series of simulations (it 
requires more than 1.2 Gbytes of memory). 

The simulations converge to a result similar to that of Klein et al. (|1999|) apart from 
the one bar fragment. In all three simulations the bar fragment appears to form when the 
binary components have approached each other enough so that they start to spiral in towards 
their eventual merger (see discussion in ^3.3.2|) . This may imply that the bar fragment is 
produced by a tidal disruption on the bar by the first binary encounter. However, SPH with 
reduced noise and a medium resolution is still not able to prevent this artificial fragment 
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from forming. 

The fact that the peak densities of the high-resolution simulations are lower than those 
of the low-resolution simulations is due to heating being applied at a lower density in the 
former. For a comparison, we have conducted a simulation with A^=160,000, A'^„=100 and 
po = 5 X 10"^'^ g cm~^. The peak density at the end of this simulation (fe„d ~1.277 tfj 
- bottom right hand panel of Fig. has reached Ppeak = 1-9 x 10^^'^ g cm~^, a value 

much higher than those of the low resolution simulations with the same amount of adiabatic 
heating, at the same time. 

In the simulation with A^=160,000, A^„=100 and po = 5 x 10~^^ g cm~^ the bar fragments 
earlier than that with A^=80,000, iV„=50 and po = 5 x lO^^^ g cm"^ ( ^3.3.3|) and it produces 9 
fragments. Some of them merge with the binary components. Again each spiral tail contains 
two lumps. Reducing the noise for the high resolution simulations with po = 5 x 10^^^ g 
cm^^ has not prevented the bar from fragmenting into quite a few fragments. We conclude 
that for these simulations convergence is achieved from the simulation with A^„=80,000 and 
A^=50 (Fig. 13.5(1 . but clearly not to the right solution as given in Klein et al. (|1999|) . 

Again, we conclude that, with higher resolution, convergence is achieved earlier, as it was 
in the isothermal simulations. More importantly, we conclude that if the Jeans condition is 
not obeyed then by having the noise decreased we cannot obtain the expected result. Only 
when the Jeans condition was obeyed did we obtain a result close to that of Klein et al. 

With medium resolution, we have shown that we achieve convergence to the expected 
solution (or close to it) with A'„=100 both in isothermal simulations and simulations with 
adiabatic heating. To obtain the exact solution in SPH simulations with adiabatic heating, 
one would need to increase the resolution as well. 

3.5 Changing the kernel 

Since in most of the above simulations with adiabatic heating the bar has fragmented, one 
could argue that the M4-kernel we have used is unable to prevent artificial clustering of 
particles and therefore artificial fragmentation (see Fig. 12. H and the discussion on the gradient 
of the kernel in ^4.3.2|) . We have experimented with another kernel that possesses most of 
the properties of the M4-kernel: it has compact support and it is truncated at r = 2h, it 
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Figure 3.12. Left : Radial density profile of an isolated particle of unit mass and smoothing length h 
using the kernel of Eqn. 13.11 Right: The gradient of the radial density profile of an isolated particle 
of unit mass and smoothing length h using the kernel of Eqn. 13.11 

has almost the same value at r = 0, its value at r = /i is one quarter of the value at r = 0, 
it contains half of the mass in the (0, h) range and the other half in the {h, 2h) range, its 
gradient has almost the same value at r = 0, there are ~ 47 neighbours within h for a uniform 
distribution of particles. This kernel, its derivative and volume integral are given by: 



W{s) 



W (s) 



0, s ^ 2, 



0, O 2, 



(3.1) 



(3.2) 



and 



W*(s) = —{ 
^ ^ 16 ^ 



(3.3) 



3s^ - 15s^ + 20s^, ^ s ^ 2; 
16, s^2. 

The fact that the second derivative of the kernel is always positive, W (s) = ensures 
that the magnitude of the gradient of the kernel decreases monotonically with increasing s 
and therefore, there is no inversion in the direction of the hydrostatic force. The fact that 
ensures that there is not going to be formed artificially any pairs of particles due 
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to a vanishing hydrostatic force. The density profile and its gradient for an isolated particle 
of unit mass and smoothing length h, using Eqns. 13. II fc IS?^ respectively, are shown in Fig. 

Em 

We have used this kernel in simulations of a rotating, spherical, uniform-density, isother- 
mal cloud with an m = 2 perturbation. In fact, we have conducted two different simulations: 
in one we have included the new kernel only for calculating the hydrodynamics (we use the 
M4 for the gravity calculations) and in the other we have used the new kernel for both gravity 
and hydrodynamics. This way, we can make a clearer comparison with the performance of 
the M4 kernel. We have used 80,000 particles initially taken from a lattice. We have included 
adiabatic heating starting at po = 10~^^ g cm~'^. The results should be compared with the 
corresponding simulation in M3.3.1I fFig. 13. 2|) . Both simulations were ended at tend ~1.295 

tff (Fig. Em. 

The binary forms at the expected time (t;,j„, ~1.25 tff). The bar produces several frag- 
ments in random positions and there are also lumps in the spiral tails (one in each). 

The peak density at the end of both simulations is Ppeak = 3.8 x 10^^^ g cm~^, almost 
twice the value of the simulation of ^3.3.11 This implies that the difference is due to the 
hydrodynamics of the new kernel. In particular, the new kernel appears to produce more 
centrally condensed objects than the M4 kernel. 

This may also explain the reason for the earlier bar fragmentation {tfrag ~1.275 tff) in 
the simulation where the new kernel was only used for calculating hydrodynamics. In the 
other simulation, the bar fragmented at the same time as the simulation of ^3.3.11 (i/rag 
~1.280 tff). 

Finally, the fact that the simulation with the new kernel used only for hydrodynamics 
produces less fragments than the simulation with the new kernel used for both gravity and 
hydrodynamics, should be attributed to the effect the new kernel has on the gravity calcu- 
lation. In particular, it seems that the new kernel emphasises density enhancements in the 
bar so that they become self-gravitating faster. 

We conclude that a kernel that theoretically should prevent artificial clustering of particles 
from occurring appears to be more problematic than the M4. This does not prove that the 
M4-kernel is the ideal kernel for Star Formation simulations. It just indicates that the artificial 
fragmentation of the bar with the M4 is not only due to it allowing artificial clustering of 
particles. Further extensive comparison studies should be made with the use of other kernels 
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Figure 3.13. Column density plots for a cloud of 80,000 particles using the new kernel (po = 10~^^ 
g cm~^). The linear size of these plots is 0.004 pc. The colour table has units of 1.18 x 10^ g cm~^. 
Top: Column density plot at the end {t = 1.293 tff) using the new kernel only for hydrodynamics. 
Bottom : Column density plot at the end {t = 1.295 tff) using the new kernel for both gravity and 
hydro dynamics . 
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due to its large size 



Figure 3.14. Column density plot for a cloud of 80,000 particles at the end with particles initially 
taken from a settled distribution (po = 5 x 10^^** g cm~^). The time is t — 1.31 t//. The linear size 
of this plot is 0.004 pc. The colour table has units of 1.18 x 10^ g cm~^. 

to identify those that operate best with self-gravitating SPH codes. 

3.6 Settled distribution 

Finally, we have conducted a simulation with 80,000 particles initially taken from a settled 
distribution. Adiabatic heating starts at po = 5 x 10~^^ g cm~'^ in order for the Jeans 
condition to be obeyed. The results should be compared with the corresponding simulation 
in ^3.3.2I (Fis. I3.3jl . Fig. 13.141 is a column density plot at the end of the simulation (at tend 
= 1-312 tff). 

The binary forms at the expected time but the bar fragments into 3 objects, formed in 
random positions. In fact, one of them forms from the merger of two smaller ones, something 
which was never observed in the simulations where particles where initially taken from a 
lattice. Having more bar fragments and them being in random positions, shows that a 
settled distribution clearly contains more noise than a distribution of particles taken from a 
lattice. This is also suggested by the fact that the peak density at the end, Ppeak = 7.5 x 
10^12 g cm^'^, is lower than that of the corresponding simulation in ^3.3.21 
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The bar fragmentation is delayed with the settled distribution {tjrag = 1-298 i//)- This 
indicates that it takes more time for a self- gravitating object to form in the bar due to 
the cloud containing more noise. In fact, the initial noise seems to provide an extra means 
of (turbulent) support to the bar. This is not necessarily a disadvantage for realistic Star 
Formation simulations, as the initial conditions for Star Formation in nature are far from 
being smooth and without noise. In chapter [31 we shall use settled distributions for our 
initial conditions. 
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Particle Splitting 

4.1 Jeans condition 

Several authors have recently discussed the significance of results based on numerical simu- 
lations that use either Eulerian or Lagrangian formulations, in particular the ability of these 
codes to prevent the non-physical growth of numerical perturbations, and their ability to 
resolve all the structure formed and therefore produce reliable and realistic results. 

Truelove et al. H1997I 11998 1) describe the need for a condition to regulate the linear 
size of their grid. They have concluded that the linear size of their smallest grid must 
always be smaller than a quarter of the local Jeans length. By setting this Jeans condition, 
they suppress the formation and propagation of artificially induced perturbations that could 
otherwise corrupt their results. 

In SPH, such a condition is also needed for similar reasons. As the SPH particles are 
allowed to move with the fluid, the properties of the fluid are kept updated. Therefore, 
the artiflcial growth of perturbations is inhibited, provided that the SPH calculations give 
accurate estimates for these properties, i.e. there is adequate sampling of the fluid with 
enough particles of similar properties within each kernel. Bate & Burkert (|1997j) have argued 
that SPH, with its current formulations, does not necessarily fulfll this condition. They 
demonstrate the need for a Jeans condition for SPH, viz. that the local Jeans mass should 
be resolved at all times. By this, they mean that there should always be enough sampling 
points in a clump near to or above the locally deflned Jeans mass. 

They also give a number of complementary rules for the formulation of the Jeans condition 
in order to obtain reliable results in fragmentation simulations. They suggest smoothing the 
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hydro dynamical forces at a scale similar to the gravity softening (i.e. e = h), as they find that 
if e < /i then artificial fragmentation is induced, while for e > h fragmentation is inhibited. 

Whitworth H1998|) has also looked into the definition of the Jeans condition for SPH and 
his findings are in very good agreement with those of Bate & Burkert (11997(1 . In his analysis, 
he introduces 3 different masses: the minimum resolvable mass by SPH, Mmim i-e. the mass 
within the radius of a kernel, with Mmin = -^n 'mptd, where Nn is the number of neighbours 
within a kernel (equal to ~50 in 3-dimensional simulations) and uiptci the mass of each SPH 
particle^; the mass of a proto-condensation, Mq; and the local Jeans mass, Mj, defined for the 
gas confined in and around Mq. If Mq ^ M^m the proto-condensation is resolved. For Mq ^ 
Mmin, the proto-condensation is unresolved. If Mq > Mj it is unstable against collapse. For 
Mo < M J, it is stable. 

He concludes that SPH is treating resolved proto-condensations (Mq S> Mmin) ade- 
quately'^. Problems may arise when they are unresolved. In this case, unresolved proto- 
condensations that are Jeans unstable (Mj < Mq <C Mmin) are not allowed to collapse as 
fragmentation is inhibited, while for unresolved proto-condensations that are Jeans stable 
(Mo <^ Mmin & Mo < Mj) artificial fragmentation is induced. He shows that artificial 
fragmentation is prevented provided that the Jeans mass is resolved, Mmin <Mj, and the 
interpolating kernel is sufficiently centrally condensed. This condition gives a very strong con- 
straint, as it does not allow any unresolved proto-condensations to form (neither stable nor 
unstable ones) and thus eliminates all possible problems that could arise in a fragmentation 
calculation. 

Bate & Burkert (|1997(1 give an even stronger Jeans condition which states that 



The proof for this is not robust, as it is based on qualitative evidence from weak convergence 
of a series of numerical simulations to a certain result. However, we will use this strong Jeans 
condition, as, even in the limiting case of 2 Mmin ~ Mj, the Whitworth (|1998|) condition is 
still satisfied. 

^It is assumed that the simulation is implemented with particles of equal mass. 

■^There is a small under-estimation of the time scale for the growth of a possible Jeans instability, which 
becomes significant only for Mo ~ Mmin- 



2Mmin < Mj. 
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Using Eqn. IA.81 the limiting case of the Jeans condition (Eqn. I4.1() becomes 



2N, 



Ntotal 



2Nnmptci = 2M„i„, < Mj 



^3^5/2 



(4.2) 



n 



1/2 



'max 



where M.totai and '^^totai are respectively the total mass and the total number of particles in 
the simulation, and Pmax is the maximum resolvable density. For a clump of given mass 
and temperature, the maximum resolvable density is a function only of the total number of 
particles 



The objective of numerical simulations of star formation is to approach stellar densities, 
i.e. to achieve as high maximum densities as possible. Eqn. 14.31 clearly sets an obstacle 
against the implementation of this objective. For example, for resolving fully a simulation 
involving an isothermal (cq = 0.17 km s~^) clump of IMq to a density of 10^^*^ g cm"^ then 
~ 1.8x10^ particles are needed, which is at the limit of present-day computer capabilities. 

An alternative would be to redirect resources only to regions of particular interest in an 
existing simulation. This can be achieved by increasing the number of particles locally in 
order to maintain the validity of Eqn. 14. II in regions approaching the resolution limit, whilst 
retaining the coarse resolution in resolved regions. This way, we can balance the need for 
higher resolution against present-day computer capabilities. 

We have invented a method to implement this. All particles in a region of interest are 
split. We give this method the obvious name "Particle Splitting". The development and 
testing of this method constitute one of the primary aims of this work. Since we can use this 
method at several levels every time the resolution limit is reached during a simulation, we 
have introduced the notion of simulations of increasingly high resolution nested inside the 
original coarse simulation. From this, we have named all simulations to which this method 
is applied "Nested Simulations". 

With particle splitting we can also address another point made by Bate &: Burkert H1997|) : 
we can follow the detailed evolution of all fragments as well as the global evolution of the 
simulation. The ability to do so is of course constrained by the time-step. 



Pmax 




(4.3) 
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4.2 Particle splitting: Concept 

Our aim is to increase the number of particles in a small sub-region of the computational 
domain of an existing simulation, just before this region reaches its resolution limit (Eqn. 14. 
This way, we will be able to continue the simulation at a higher resolution, but only where this 
is really necessary. The shape of the sub-region will depend on the geometry of the problem 
(e.g. cylindrical in simulations involving flattened structures), but the simulation will always 
be fully 3-dimensional. In the sequel we shall refer to particles in the high resolution region 
as fine particles, and particles in all other regions as coarse particles. We shall also refer to 
a simulation as a fine simulation, if it includes fine particles, and as a coarse simulation, if it 
includes no fine particles. 

The method will be applied at time tgpi, when significant - but always resolved - structure 
has formed in the coarse simulation, just before this structure reaches its resolution limit. 
We will stop the coarse simulation and decide the position, shape and size of the sub-region 
manually. This involves choosing the appropriate co-ordinates for the sub-region, so that it 
contains all the significant structure. The initial conditions for the fine simulation will be 
interpolated from the coarse simulation. At tspi, each coarse particle in the sub-region will 
be replaced by 13 new equal-mass fine particles distributed symmetrically round the coarse 
particle's position f ^4. 3. 1(1 . These fine particles will be given velocities interpolated from the 
velocities of the coarse particles they replace f ^4.3.5|) . Using the appropriate subroutines of 
the code (chapter they will also be given values for all the physical quantities involved in 
the SPH equations, such as density f ^4.3.2l fc l4.3.3() . smoothing length ( ^4.3.4|) . temperature, 
acceleration, time-step information ( ^4.3.6|) . After this initialisation, the fine simulation, 
with only fine particles inside the sub-region and only coarse particles outside the sub-region, 
will start. In subsequent time-steps, if a coarse particle is found to cross the sub-region's 
boundary, then it will be split on-the-fiy into 13 fine particles and with a procedure similar 
to the initialising one, these fine particles will be given velocities, densities, temperatures, 
smoothing lengths, accelerations, time-steps, and they will immediately become active. On 
the other hand, if a fine particle exits the nested sub-region, it will continue being active. The 
positions, velocities and accelerations of the coarse particles outside the sub-region will evolve 
together with those of the fine particles. In effect, the outside coarse particles will provide 
the boundary conditions to the fine simulation. This means that the boundary conditions 
will be as exact as those of the coarse simulation. 
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A possible problem that we will have to deal with is the interaction at the boundary of 
two different populations of particles: the massive and extended coarse particles, and the 
light and compact fine particles. This can cause interpenetration and mixing of the two 
populations and an artificial blurring of the boundary. We solve this problem by adjusting h 
so that the kernel contains a fixed multiple of the mass of the central particle f ^4.3.4|l . rather 
than a fixed number of neighbours. 

An alternative implementation of particle splitting involves setting a threshold density 
above which particles are automatically split. The advantage of the latter method is that we 
don't have to stop the simulation to decide the co-ordinates of the sub-region. All splitting 
happens on-the-fiy. We call this version of the new method "on-the-fly splitting" as opposed 
to the former version which we call "nested splitting" . 

On-the-fiy splitting is our preferred version of the new method. However, we will discuss 
both methods. We will first define the positions of the fine particles with respect to the posi- 
tions of the coarse particles they replace. This will involve fine-tuning each coarse particle's 
h and the distance between the 13 fine particles. Therefore, this will require adjusting the 
density profile of the configuration of the 13 fine particles to the density profile of the coarse 
particle they replace as well as trying to match the density profile of coarse particle distribu- 
tions with the density profile after particle splitting. We will then illustrate the difficulty of 
simulating the boundary between fine and coarse particles. Subsequently, we will describe the 
new method for finding h for all particles both fine and coarse. The next section will finish 
by describing the way we assign other physical and numerical (e.g. velocity, acceleration, 
temperature and time step information) properties to all fine particles. 

We will then perform some tests (^QJ in order to validate the performance of the method. 
We will simulate uniform collapse, where we will also use sink particles ( ^4.4.1|) . We will show 
the efficiency of the new method for calculating h by applying it to the simulation of a stable 
isothermal sphere ( ij4.4.2|l . Finally, we will apply particle splitting to a collapse simulation 
of a rotating, uniform-density, isothermal cloud with an m=2 perturbation (see chapter ^ 
to utilise the new method in a more realistic application f ^4.4.3|) . showing its efficiency in 
reproducing accurate results, as well as economy in terms of computational cost. 

In the next chapter, we will apply particle splitting to clump-clump collision simulations. 
Firstly, we will extend previous simulations of high mass cloud collisions to quantify the 
benefits of the new method. Secondly we will investigate a new part of the parameter space. 
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0.0 
Van 



Figure 4.1. Graphical representation of the configuration of 13 fine particles, in two dimensional 
projection (left) and three dimensions (right). 

looking for realistic fragmentation mechanisms in collisions between low-mass clumps. 



4.3 Particle splitting: Implementation 
4.3.1 Positions for the fine particles 

We replace a coarse particle with a configuration of 13 particles as shown in Fig. 14.11 A fine 
particle is put exactly at the position of the coarse particle, while 6 particles are put on the 
vertices of a hexagon centred on the position of the coarse particle. The remaining 6 particles 
are put on the vertices of two equilateral triangles parallel to the plane of the hexagon but 
on either side. All particles are put at equal distances rj from their nearest neighbours^. 
Such a lattice is the simplest possible arrangement that has minimum interstitial volume 
()Kittel 1962|) . In fact, it is the primitive cell of a face-centred cubic (cubic closed-packed) 
structure (see Fig. 14.11 left panel). Specifically, Table H?T] gives the co-ordinates of the 13 fine 
particles at unit distance away from each other, in the reference frame of the coarse particle 
they replace. 

The value of rj is of great importance as well as the value of the smoothing length of 
^We will deal with the value of rt shortly. 
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Table 4.1. Co-ordinates of the 13 fine particles at unit distance away from each other, in the reference 
frame of the coarse particle they replace. 

the fine particles, hi. They should be fine-tuned to minimise the deviation between the 
configuration of 13 fine particles and their parent coarse particle. Apart from reproducing 
the density profile of the parent coarse particle, the 13 fine particles should be positioned in 
such a way that they experience the same accelerations as if their parent coarse particle was 
present. Finally, we should take care that any density fluctuations input to the fine region 
after particle splitting are kept to a minimum. We will now calculate the density profile of 
the configuration of 13 fine particles in isolation as well as estimate any density fluctuations 
in a particle distribution due to the application of particle splitting. 

4.3.2 Density profile of the configuration of 13 fine particles 

In order to calculate the mean density profile of such a configuration of particles, we consider 
a sphere of radius r centred on the central fine particle. We define the distance d from the 
position of any one of the 12 outside fine particles to any point on this sphere as a function 
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The figure is provided separately 
due to its large size 



Figure 4.2. Two dimensional projection of the model used for the calculation of the mass distribution. 
The two concentric circles represent the one and two smoothing length spheres of an outside fine 
particle. The sphere with radius r is centred on the central fine particle which is at a distance r, away 
from each outside fine particle. 
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of r and the angle 9 (see Fig. 14.2(1 as 

d = (r^ + rf - 2rriC0s{e)) ^ . (4.4) 
The mean density on the surface of the sphere is then evaluated using 



p{r)4:TTr'^5r = ^"^ j rmKr^WKu [^jl^ 5r dA 



p{r) = Qrriihf I Wma ( ^ ] sin{e)de, (4.5) 



hi 



where the first equation gives the total mass swept up by a shell of radius r and infinitesimal 
thickness 5r. We have used A = 2'irr'^ sin{0) , the surface area of the sphere that is limited by 
the circle defined by angle 9 (produced when d, the dashed line on Fig. 14.21 is rotated about 
Tj through an angle (p = 2tt). rrii is the mass of one of the fine particles (rrij = M/13 where 
M is the mass of the coarse particle) and /ij its smoothing length (/ij = (i//13)^/'^, where H 
is the smoothing length of the coarse particle). We have not taken into account the mass of 
the central fine particle yet. Since we are averaging over all angles we have multiplied by 12 
in order to account for all the outside fine particles. 
The integral of Eqn. 14.51 becomes 



p(r) = Gruih- / Wm4 

f-^max 



dp, (4.6) 



hi 

where we have substituted p = cos{9). We have used the minus sign of the dcos{9) term to 
exchange the limits of the integral. 
If we write Wm4:{s) in the form 

Wm4:{s) =wo + Wis + 1028"^ H h WpS^ + . . . 

then Eqn. 14.61 takes its final form 

-pir) = Qm^hf r' r r^ + rf -2rnp Y/'] 

J ^J.=^lm.ax I p—Q V ""i / I 

with the sum being the polynomial form of the M4- kernel (Eqn. 12.5(1 . 

The value of pmax has not been defined yet. There are four different cases that we must 
take into account for the value of Pmax, depending on the value of r. In particular: 
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Figure 4.3. Linear density profiles of a coarse particle in isolation (solid curve) and of the ensemble 
of 13 fine particles that replace it (dashed curve). The distance of the 12 outside particles from the 
centre is set to be exactly 2hi in the left panel and 1.9hi in the right panel, and it is indicated by the 
vertical lines. 

1. If the sphere centred on the central fine particle intercepts the two smoothing length 
sphere of an outside fine particle but not the one smoothing length sphere {hi ^ r| ^ 



(Fig. 14. 2() . Therefore, in this case. 



and for Eqn. 14. 71 we use the second part of the M4- kernel (Eqn. 12.5(1 for 1 ^ s ^ 2 and 

2. If the central fine particle lies inside the 2h spheres of the outside fine particles and 
outside of the h spheres without intercepting the h spheres {hi ^ r, — r and rj + r ^ 2/ij, 
Fig. Ol), then 

^max — and for Eqn. 14.71 we use the second part of the M4-kernel 
(Eqn. 12. 5|) for 1 ^ s ^ 2 and ^imax = — 1- 

3. If the sphere centred on the central fine particle intercepts both spheres {\ri — r| ^ hi, 
Fig. 14. 2|) . then the integral of Eqn. 14.71 breaks into two parts 



2/ii, Fig. 1321), tlien 9, 



max — 



^2, where 
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p{r) = Qmih\ 



-3 




where we use 



(Fig. 1121) so that 



— ^ — J Y'^ 



7i = cos 



2rrj 



+ rf - hf 



2rrj 

and fimax = /U2- For Eqn. 14.81 we take Wp and p from the first part of the M4-kernel 
(Eqn. 12.51 for ^ s ^ 1) for the first sum, and from the second part of the M4-kernel 
(Eqn. 12.51 for 1 ^ s ^ 2) for the second sum. 

4. If the central fine particle hes inside the 2h spheres of the outside fine particles and 
outside of the h spheres but it now intercepts the h spheres (hi ^ — r and n + r ^ 2hi, 
Fig. 14. 2|) then we use Eqn. 14.81 with 6max = t^-, thus pmax = — 1- 

The analytical calculation of these two integrals has now become trivial as, in every case 
we have to deal with integration of polynomial functions. In all other cases, the integral 
becomes equal to due to the M4- kernel having compact support. 

There are two parameters of the problem which we have not dealt with yet. Namely, the 
density of the central fine particle and the value of rj, the distance from the central to the 
outside fine particles. The calculation of the former is based on Eqn. 12.71 with = the 
position of the central fine particle, = M/13 the mass and hi = (/7/13)^/^ the smoothing 
length of the fine particles. The total density p{r) for the ensemble of 13 fine particles is 
derived when we add the density of the central particle to the result of Eqns. 14.71 and/or 
14.81 This is illustrated with the dashed curve on the left panel of Fig. 14.31 with the outside 
particles at a distance of rj = 2/ij. 

Vi is a free parameter in this problem and its value can be determined by minimising the 
integral 

ATir'^\p{r) - p{r)\dr, (4.9) 

r=0 
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Figure 4.4. Gradient of the density profiles of a coarse particle in isolation (solid curve) and of the 
ensemble of 13 fine particles that replace it (dashed curve). The distance of the 12 outside particles 
from the centre is set to be exactly 2/ij in the left panel and 1.9/ii in the right panel, and it is indicated 
by the vertical lines. 

where p{r) is the density of the coarse particle"^ and p{r) is the density of the configuration 
of 13 fine particles that replace it. This integral gives the volume average of the absolute 
difference in the density estimates given by p{r) and p{r) as defined above. It is calculated 
numerically using the extended trapezoidal rule 



where fi are the values of the function / evaluated at Xj (i = 1, 2, . . . , N) and h = Xi — Xi-i 
{i = 2,3, . . . , N) the constant step (jPress et al. 1990|1 . 

The minimum of integral 14.91 gives = 1.9^j (i.e. 95% of 2/ij). The right panel of Fig. 
14.31 shows p(r) and p{r) (solid and dashed curves respectively) when rj = 1.9/ij. The left 
panel of Fig. 14.51 shows a zoom on this, for r between and 0.2. 

Examination of Fig. 14.31 shows that the shape of the density profile for the ensemble of 
fine particles is significantly different compared to the density profile of a coarse particle. 
The dip that appears in the former is due to the geometry of the fine particle configuration, 
i.e. the density drops with distance as we move away from the central particle before it rises 
■'its calculation is similar to that of the central fine particle. 
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Figure 4.5. Left : Zooming on the right panel of Fig. 14.31 Density profiles of a coarse particle in 
isolation (solid curve) and of the ensemble of 13 fine particles that replace it (dashed curve). The 
distance of the 12 outside particles from the centre is set to be 1.9hi. Right: Gradient of the density 
profiles of a coarse particle in isolation (solid curve) and of the ensemble of 13 fine particles that 
replace it (dashed curve). The distance of the 12 particles from the centre is set to be I.Shi, and it is 
indicated by the vertical line. 

again due to the outside particles. 

We have also checked the behaviour of the density gradient. Fig. 14.41 shows the density 
gradient of a coarse particle (solid curve) and of the group of fine particles that replace it 
(dashed curve), calculated numerically^ for rj = 2/ij and rj = 1.9/ii (on the left and right 
panel respectively) . One non-physical property of the M4- kernel is that the density gradient 
vanishes as r — > 0, and therefore if two particles are very close the repulsive hydrostatic 
force between them is very small. This means that there will be a length scale below which 
artificial clustering of particles may occur. We would expect the density gradient of the 13 fine 
particles that will replace a coarse particle to behave at least similarly to that of their parent 
coarse particle. However, the density gradient of the group of fine particles has two minima 
(dashed lines of Fig. 14. 4 j) which is worrying. Firstly, this means that there are going to be two 
length scales between which particles will be pulled towards the centre. Secondly, the smaller 

^ [^] i ~ (Pi~Pi-i)/'*5 where pi is the value of the density evaluated at a;^ (i = 1, 2, . . . , A'") and h = Xi—Xi-i 
{i = 2,3, . . . , N) the constant step. 
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Figure 4.6. Left : j immediately after the splitting vs. Vi. Right: the density distribution of the 
initial settled box (before splitting), evaluated on the particle positions. 

of these two length scales is a point of stable equilibrium, i.e. artificially clustered groups of 
particles may form and persist. Fortunately, it seems that there is an intrinsic geometrical 
constraint in the process of replacing coarse particles by fine ones. The fine particles are 
always positioned outside the outer of the two length scales where the hydrodynamical forces 
are reversed in direction. This means that our fine particles will always be positioned at 
a region where they will experience outward hydrodynamical forces. Problems may appear 
only if they are perturbed inwards and end up in the non-physical region. It appears that 
this does not happen in practice. 

In order to explore if a group of fine particles can produce a density gradient that behaves 
similarly to that of a coarse particle, we have conducted another short parameter search for 
the best value of rj. For rj ^ 1.3/ij (right panel of Fig. 14.51 for rj = 1.3/ij) the density 
gradient appears to take a shape similar to that of a coarse particle. However, we found that 
for = 1.3/ii, integral 1)4.9(1 gives a value for the absolute difference in the density estimates 
twice as large as the = 1.9/ij does. 

4.3.3 Density stability with particle splitting 

In this section we demonstrate that particle splitting does not significantly affect an existing 
density distribution of coarse particles. We have produced a uniform distribution of 500 
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Figure 4.7. Density distribution of the box after splitting for = 1.5hi (left panel) and = 1.9hi 
(right panel), evaluated on the particle positions. 

particles in a box using periodic boundary conditions. We have evaluated the density both 
on the particle positions and on a square grid. Using the method described in ^4.3.11 we then 
replaced each particle with 13 fine ones randomly oriented and we tried to reproduce the 
same uniform density allowing for a short settling of the fine particles. We have discovered 
that the value of rj, the distance of the 12 outside fine particle from the central fine particle, 
can influence our results significantly. 

In particular, the amount of settling required after particle splitting for the density dis- 
tribution to be identical to the initial settled box depended on this distance r,. We have used 
the ratio of the variance over the mean of the density distribution, j, as the quantity that 
determines if our distribution is settled or not. We have settled until ^ < 0.01. The left panel 
of Fig. 14.61 shows the dependence of j on immediately after the splitting. The right panel 
illustrates the density distribution of the initial settled box (before splitting). The minimum 
value of J is for rj between 1.5/ij and 1.6/ij, with 1.5hi having the absolute minimum in both 
variance and mean of the density. Fig. 14.71 shows the density distribution of the box after 
splitting for rj = 1.5/ij (left panel) and rj = 1.9/ij (right panel). Thus, for rj = 1.5/ij the least 
settling is required. This applies for both the evaluation of density on particle positions and 
on the square grid. 

Therefore, we finally decided that rj should neither be equal to 1.9/ij nor 1.3/ij (c/. ^4.3.2|1 
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Figure 4.8. Density profile (left) and density gradient (right) of the ensemble of 13 fine particles 
(dashed lines) that replaced a coarse particle (solid lines). The outside fine particles are at — 1.5/ii 
from the centre. The vertical lines indicate the position of r^. 

but 1.5/ij. This is because it is the distribution of an ensemble of particles that we are 
particularly interested in, and not the density profile of a single particle. Any over-evaluation 
of the density around the position of the coarse particle (c/. Fig. 14.81 left panel) is expected 
to be treated appropriately by the hydrodynamical forces (dispersal of the ensemble of fine 
particles that replace it). The shape of the density gradient reassures us that this will happen 
(right panel of Fig. I4.8jl . 

In ^4.4.31 we will discuss the effects of particle splitting on a more realistic particle distri- 
bution. 

4.3.4 Smoothing lengths for the fine particles 

As shown in Fig. 14.61 & 14.71 there is an increase in the mean density of a uniform distribution 
of 500 particles after particle splitting is applied, even for the case of rj = 1.5/ij when this 
increase is minimum (i.e. minimum settling is required for the fine distribution to obtain 
again unit mean density). In particular, for this case the mean density is increased by 5.04%. 
Furthermore, the calculation of the number of neighbours for each fine particle does not 
always give the right result (~ 50). 

In order to explore this further we have conducted another similar test. We have simulated 
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a stable isothermal sphere in isolation and have applied particle splitting (in particular the 
nested splitting version) to its central region (the initial conditions are given in H4.4.2I and 
shown in Fig. I4.14|) . We have found that it is not only the number of neighbours that 
gives wrong results, but also the density inside and outside the fine sub-region is incorrectly 
modelled after splitting is applied. The reason for this is that there now exist particles of 
different mass in close contact. In particular, the method for finding the smoothing length 
for all particles is based on counting the number of their neighbours ( ^2.6|) . This makes 
fine particles just inside the sub-region boundary, which have coarse particles as neighbours, 
look for some of their neighbours in a region of lower resolution. Therefore, their smoothing 
length is over-estimated and their density under-estimated. The coarse particles just outside 
the sub-region boundary have fine particles as neighbours and look for their neighbours in a 
region of higher resolution. Thus, their smoothing length is under-estimated and their density 
over-estimated (top panel of Fig. 14.91 cf. with top panel of Fig. I4.14j) . However, there is an 
even more severe boundary effect disturbing the evolution of the fine simulation. Since the 
sub-region boundary is fixed in space, the fine particles should always be inside the boundary. 
Nevertheless, with the implementation of nested splitting fine particles penetrate through the 
non-moving coarse particles and end up at the other side of the boundary (bottom panel of 
Fig. 14.91 cf. with bottom panel of Fig. I4.14() . The gravitational field around the boundary 
fine particles is not balanced as they are in contact with more massive particles from one 
side. This make them move to the other side of the boundary. In particular, they move 
through the gaps between the boundary coarse particles, occupying the potential energy 
gaps between those particles. This creates an expanding shell of fine particles similar to 
an artificial rarefraction wave which will eventually corrupt the evolution of the isothermal 
sphere even if the sphere may finally evolve to another stable equilibrium state. 

Therefore, the new method clearly needs a special feature to eliminate these boundary 
effects. A possible solution would be to define a region on either side of the boundary, where 
both fine and coarse particles would have their smoothing length evolution constrained in 
order to create a smooth transition from the low to the high resolution and vice versa. 
However, this idea would be rather complicated to implement as the position and velocity of 
particles with respect to the boundary would need to be calculated at every time-step. 

We have found a simpler method. It involves calculating the smoothing length of a particle 
by specifying the total mass of neighbours, rather than the number of neighbours. Specifically, 
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the smoothing length of each particle is such that it contains ~50 times its mass. We think 
this is more appropriate in our case where there exists mixing of particles with different mass. 
Away from the boundary, finding h is exactly the same as before: there are ~ 50 equal mass 
neighbours. Close to the boundary, though, a coarse particle appears as 13 fine ones when h 
is calculated for a fine particle, while it takes 13 fine particles to add one effective neighbour 
to a coarse particle. This method is valid as the volume of a coarse particle (oc H"^) is not 
significantly altered when the coarse particle is split with rj = 1.5/ii ( N4.3.2j) ^. Close to the 
boundary, coarse particles have more than 50 neighbours and fine particles have less than 
50, but both their smoothing lengths and densities are correctly modelled. In its numerical 
details, the method is implemented exactly the same way as before ( ^2.H|) . but the variable 
that we now use as a criterion for the acceptance of a trial smoothing length is the enclosed 
mass and not the number of enclosed neighbours. The notion of a particle always tracing 
constant mass is maintained and SPH with different mass particles preserves its Lagrangian 
character. The test simulation for the evolution of a stable isothermal sphere is repeated 
f M.4.2|) and shows that the new method for calculating h greatly improves the treatment of 
the boundary as well as preventing the growth of numerical perturbations that may disturb 
the evolution of the stable isothermal sphere. 

4.3.5 Velocities for the fine particles 

To evaluate the velocity of each fine particle we use Eqn. 12.61 

v(r,) = Y.^.^Hr^WM4 (^^) > (4-10) 

where Vj, rj and Hi are the velocity, position and smoothing length, respectively, of coarse 
particle i, and rj is the position of fine particle j. We do not sum contributions from coarse 
particles which are more than 2Hi away from fine particle j (i.e. '"^^j-'"' > 2). Since each fine 
particle is contained within the 2H radius of its parent coarse particle, the sum of Eqn. 14.101 
is over the ~50 neighbours of its parent particle. Obviously, the largest contribution comes 
from this parent particle, but Eqn. 14. 101 guarantees that the velocity field around each coarse 
particle is accurately passed to the fine particles that replace it. 

In all the tests to which particle splitting has been applied ( ^4.4|) . the velocity field is 
accurately represented by the fine particles after splitting occurs. The new method also 
^For the first trial value of hi = the configuration of 13 fine particles has a radius of ~ 1.06-ff. 
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Figure 4.9. Evolution of a stable isothermal sphere after t = 9.07 tff, when nested splitting was 
applied within a radius of 3 x 10^^ pc, without implementation of the new method for calculating h. 
The green points show fine particles and the black points coarse particles. Top: Radial density profile 
of the isothermal sphere. The red line indicates the solution of Eqn. 12.541 Bottom : Thin equatorial 
slice (Az = 4 X 10^'^ pc) of the isothermal sphere showing the boundary between the coarse and the 
fine particles. 



90 



CHAPTER 4. 



PARTICLE SPLITTING 



Starting step n 
t_now = t(n) 




Second half step 



Advance system to 
t_now 



Find active particles 



Calculate acceleration 




As in first half step 



MPT: find time bins 



Make Tree 

& 

Set h for the cells 



Calculate h for 
active particles 



Calculate density 
& 

Hydro acceleration 



T 



Calculate gravity 



Figure 4.10. The code 'step by step': Flow-chart of the algorithm that dictates the evolution of the 
fluid in time, after the on-the-fly splitting subroutines have been added. It represents the cycle, n, 
of the integration scheme, when the system advances with the time-step in the minimum time bin, 
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conserves global linear and angular momenta. 

4.3.6 Updating other fluid properties & numerical parameters 

In "nested splitting" , the exact steps we follow in implementing particle splitting are: 

• We decide on the dimensions of the fine region. We count all the coarse particles that 
are inside this region at time tgpi when the method is initiated. 

• The fine particles that are produced are put in position by randomly rotating the co- 
ordinates given in Table ETTl around the y- and z-axes. 

• The velocities of the fine particles are calculated using Eqn. 14.101 Their mass is 13 
times smaller than that of a coarse particle. Their h is calculated with the new method 
described in the previous section ( M4.3.4|) with a first trial value of hi = (i//13)^/'^. 

• The fine particles are given the same temperature as their parent coarse particle. Tak- 
ing into account that splitting is necessary in order for a simulation to be properly re- 
solved up to the density where adiabatic heating initiates (chapter [2J , this assumption 
of isothermality is valid for the fine particles, as all other particles remain isothermal 
as well. If particle splitting were necessary at later stages in regions with tempera- 
ture gradients, temperatures could be calculated by interpolation in the same way as 
velocities. 

• All fine particles are put in the minimum time bin to make sure that the method 
sustains its accuracy. They are given the minimum time-step. If the multiple time-step 
method decides that they can be moved to a higher time bin, then this will happen 
as soon as it is allowed f H2.9|) . The extra computational cost of keeping some fine 
particles in the minimum time bin ~ even though this may not really be needed - 
is minimal. Nonetheless, our decision for the initial choice of time-step for the fine 
particles is dictated by the need to obtain the highest possible accuracy. Moving all 
the fine particles and updating their properties during the first time-steps after their 
introduction to the simulation also gives the code a cushion for correcting any values 
that do not match the local fiuid properties around the fine particles. 

• The fine particles are given a splitting identifier indicating the level of splitting they 
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represent, as the method may be apphed several times in a simulation creating a series 
of nested simulations each at a higher resolution. 

• The total number of particles is updated to contain the fine particles. For book keeping 
purposes only, the central fine particle replaces its parent coarse particle in all variable 
arrays, while all other fine particles are put at the end of the arrays. 

• Finally, after the density of all fine particles is calculated, we can calculate the acceler- 
ation for these particles. For these calculations we use the usual SPH routines ( H2.7|) . 
The fine particles are now ready to take part in the next cycle of the integration scheme 
and then the fine simulation starts. 

• At the end of all subsequent time-steps, the code checks for any coarse particles that 
may have crossed the sub-region boundary during the previous time-step. Then these 
coarse particles are split on-the-fly, without stopping the fine simulation. The code just 
repeats all the above steps for the coarse particles that have just entered the fine region. 
Fig. I4.1UI shows the flow-chart of the algorithm that dictates the evolution of the fluid 
in time, after the on-the-fly splitting subroutines have been added. 

In "on-the-fly splitting" we don't have to stop the coarse simulation to apply particle 
splitting, but instead we use an automated test in the particle splitting subroutine. We 
calculate initially the density threshold, Pmax, above which the simulation stops resolving the 
Jeans mass. Pmax is given by Eqn. 14.31 Then we start the coarse simulation and particle 
splitting automatically initiates on-the-fly when particles exceed this density. 

This way, we don't have to apply splitting to particles unnecessarily, just because they 
are in spatial proximity to developing proto-condensations. Therefore, we can achieve even 
greater economy in computational cost, that can make the new method even more attractive 
and efficient. This is our preferred version of particle splitting and this is the method we 
apply to clump-clump collisions in the next chapter. In the following section we describe the 
application of both versions of particle splitting to standard test simulations. 

4.4 Tests 

We apply both versions of particle splitting to the central region of a collapse simulation, 
to see if there is propagation of boundary effects. We do the same for the central region of 
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a simulation of a stable isothermal sphere, to show that the boundary effects mentioned in 
^4.3.41 have been eliminated by the new calculation of smoothing lengths. Finally, we apply 
both versions of particle splitting to the collapse simulation of a rotating, uniform-density, 
isothermal cloud with an m=2 perturbation (c/. chapter We demonstrate that SPH with 
particle splitting gives very good results to this problem. We show that the Jeans condition 
for fragmentation provides a very strong test for the significance of numerical results. We also 
show the efficiency of the new method in terms of computational economy and in particular 
the superiority of "on-the-fly splitting" . 

In quantifying the results of our tests, we need to be able to associate these results only 
with the performance of our numerical code. In order to decrease the numerical noise input by 
the initial distribution of particles, we perform all tests using clouds whose particles are taken 
initially to be on a lattice. In the sequel, we will refer to the particles of such simulations 
as "particles initially taken on a lattice". To verify that the results of tests with particles 
initially taken on a lattice are not biased due to some preferred orientation of the initial 
lattice, we will also perform one simulation for each test where particles are taken initially 
from a "settled" distribution. Such a distribution of particles is produced when the particles 
are taken in random positions and then they are relaxed to uniform density, using the SPH 
formulation described in Whitworth et al. (|1995|1 . 

4.4.1 Collapse Simulation 

We test both versions of particle splitting at the centre of a collapse simulation. This way we 
test the method under conditions of homogeneous inflow. Both implementations show good 
results. There is very little dispersion of fine particles out of the fine region and there is a 
well-defined boundary between the fine and the coarse regions. 

The initial conditions consist of a spherical cloud of mass M = 1 Mq and radius R = 
0.016 pc. The cloud has uniform density p = 3.74 x 10"^*^ g cm~'^ and uniform temperature 
T = 7.9 K. The ratio of thermal to gravitational energy is a = 0.26 and the cloud has a 
Jeans mass of 0.2 Mq, so that it is unstable and it collapses. There is no rotation. We 
have used 10,185 particles to simulate the cloud, initially taken on a lattice. For this test we 
use our complete self-gravitating SPH code (chapter [21) . We also include adiabatic heating 
to slow down collapse. Adiabatic heating initiates at pq = 10^^^ g cm^^. It is included in 
order to prevent small time-steps from occurring. Adiabatic heating is necessary at the very 
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centre of the simulation, where matter accretes to form a central object. Outside the centre, 
isothermal collapse occurs. In this region, the collapsing gas evolves with a self-similar form, 
p oc r~^, as predicted by Bodenheimer & Sweigart H1968|) . 

Nevertheless, the time-step does get very small after some time. To prevent this, we have 
repeated the simulations using a sink particle to simulate the core of the accreted central 
object. With this, we allow accretion of more matter to the central object and we finally end 
up at a state where most coarse particles have entered the fine region and have been split. 

A sink particle is a particle that accretes all matter that enters its radius. In this case 
the sink particles have a radius of 2 x 10~^ pc. There is only one sink particle at the centre 
of each simulation. The particles that lie inside this radius stop being active (i.e. their 
properties are not followed any further). The accreted particles have their mass decreased 
to zero and their mass is added to the mass of the sink. Therefore, the region just outside 
the sink is not very well evolved, as the active particles just outside the sink radius do not 
"see" the accreted particles inside the sink radius. As a result, these active particles look 
for neighbours only outside the sink. Their smoothing length is over-estimated and their 
densities under-estimated. 

We stop all simulations when the time-step becomes less than 2 x 10~^ tjj, so that it 
is computationally inefficient to continue (i.e. in order for time to progress by 10~^ t// we 
would need 5 times the run-time up to that point). 

We present five different simulations: collapse using nested splitting without a sink, col- 
lapse using nested splitting with a sink, collapse using on-the-fly splitting without a sink, 
collapse using on-the-fly splitting with a sink, and collapse using nested splitting with parti- 
cles initially taken from a settled distribution. 

4.4.1.1. Isothermal collapse using nested splitting without a sink 

The top panel of Fig. 14. Ill shows the radial density profile of the sphere when nested splitting 
was applied with the fine region having a radius of 2 x 10"'^ pc. The green points show fine 
particles and the black points coarse particles. The red line indicates the p cx profile. 
The top panel of Fig. 14. Ill shows the end of the simulation at time t = 1.01 tjf. For about 
2 orders of magnitude in radius, we find that density complies well to the predicted profile, 
p oc r~^. Inside 3 x 10^^ pc the density increases rapidly due to the formation of a central 
core. The cloud has contracted to a radius of 6.3 x 10^^ pc. The maximum density at the 
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Figure 4.11. Radial density profile of a collapse simulation when nested splitting was applied within 
a radius of 2 x 10~^ pc. The green points show fine particles and the black points coarse particles. 
The red line indicates the p oc profile. Top: Radial density profile after t = 1.01 tff for a collapse 
simulation without a sink. Bottom : Radial density profile after t = 1.05 tff for a collapse simulation 
with a sink of radius 2 x 10~^ pc simulating the accreted central object. 
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Figure 4.12. Radial density profile of a collapse simulation when on-the-fly splitting was applied. 
The density threshold is taken to be Pmax = 8.55 x 10~^^ g cm~^. The green points show fine particles 
and the black points coarse particles. The red line indicates the p oc profile. Top: Radial density 
profile after f = 1.01 1// for a collapse simulation without a sink. Bottom : Radial density profile after 
t = 1.05 tff for a collapse simulation with a sink of radius 2 x 10~* pc simulating the accreted central 
object. 
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centre of the cloud is Ppeak = 4.8 x 10~^^ g cm~^. There are 69,861 particles in total, with 
4,973 coarse particles having been split. Note that there are two boundary effects. Some fine 
particles have their density over-estimated immediately after their parent coarse particle is 
split (i.e. their first smoothing radius is smaller than it should be). Some other particles are 
attracted by the heavier coarse particles and temporarily exit the fine region. Both these 
effects are transient, as particle identification has proven that, although the effects are static 
with time, they are produced by different (new) fine particles at each time-step. This test 
shows that the boundary effects due to nested splitting are not significant as they do not alter 
the predicted density profile. SPH with the new method for calculating h f ^4.3.4l) eventually 
remove any fluctuations input by the application of nested splitting. 

4.4.1.2. Isothermal collapse using nested splitting with a sink 

The bottom panel of Fig. 14.111 shows the radial density profile of the sphere when nested 
splitting was applied with the same radius for the fine region, but now with a sink particle 
having radius 2.0 x 10~^ pc, simulating the central accreted object. It is the end of the 
simulation at time t = 1.05 tjf. The radius of the sphere is 5.4 x 10~^ pc. Accretion to 
the centre of the cloud is not modelled properly due to the sink. The time-step does not 
decrease so rapidly and at the end of the simulation more matter has entered the fine region. 
In particular, there are 99,837 particles in total, with 7,471 coarse particles having been split. 
The radial density profile obeys the p cx profile for one order of magnitude. A small 
density under-estimation is evident outside the sink radius due to the effects described at 
the end of ^4.4.11 This test extends the conclusion that the boundary effects due to nested 
splitting are not significant in a case where most coarse particles (~ 75%) are split. 

4.4.1.3. Isothermal collapse using on-the-fly splitting without a sink 

The top panel of Fig. 14.121 shows the radial density profile of the sphere when on-the-fly 
splitting was applied. The density threshold is taken to be pmax = 8.55 x 10~^^ g cm~^. 
The green points show fine particles and the black points coarse particles. The red line 
indicates the p cx profile. The top panel of Fig. 14.121 shows the end of the simulation at 
time t = 1.01 tff. For about 2 orders of magnitude in radius, we find that density complies 
well to the predicted profile, p cx r~^. Inside 3 x 10~^ pc the density increases rapidly 
due to the formation of the central core. The density profile compares very well with the 
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density profile of the nested splitting simulation (top panel of Fig. I4.11|) . The fine region has 
smaller radius than in the nested splitting simulation, but splitting happens again within the 
isothermal region (i.e. along the p oc line), so that comparison between the two methods 
is legitimate. Again, the cloud has contracted to a radius of 6.3 x 10~^ pc, and the maximum 
density at the centre of the cloud is Ppeak = 4.8 x 10~^^ g cm~'^. In the on-the-fly splitting 
simulation there are only 46,917 particles in total, with 3,061 coarse particles having been 
split. Therefore, on-the-fly splitting, involving ~20,000 particles less than nested splitting, 
advanced to the same state as nested splitting a lot faster, requiring less computational effort 
(the required memory and the size of the output files are proportional to the number of 
particles). Note that the two boundary effects are still present. Some fine particles have 
their density over-estimated immediately after their parent coarse particle is split (i.e. their 
first smoothing radius is smaller than it should be). Some other particles are attracted by 
the heavier coarse particles and temporarily exit the fine region moving to a region of lower 
resolution. Both these effects are again transient, as particle identification has proven that, 
although the effects are static with time, they are produced by different (new) fine particles 
at each time-step. This test shows that the boundary effects due to on-the-fly splitting are 
not significant as they do not alter the predicted density profile. SPH with the new method 
for calculating h ( ^4.3.4() eventually removes any fluctuations input by the application of 
on-the-fly splitting. Comparison with nested splitting demonstrates that on-the-fly splitting 
produces the same results as nested splitting, much more efficiently. 

4.4.1.4. Isothermal collapse using on-the-fly splitting with a sink 

The bottom panel of Fig. 14.121 shows the radial density profile of the sphere when on-the-fiy 
splitting was applied with the same density threshold, but now with a sink particle having 
radius 2.0 x 10~^ pc, simulating the central accreted object. It is the end of the simulation 
at time t = 1.05 tfj. The radius of the sphere is 5.4 x 10~^ pc. Accretion to the centre of the 
cloud is not modelled properly due to the sink. The time-step does not decrease so rapidly 
and at the end of the simulation more matter has entered the fine region. In particular, 
there are 86,817 particles in total, with 6,386 coarse particles having been split. The radial 
density profile obeys the p oc r^^ profile for one order of magnitude. A small density under- 
estimation is evident outside the sink radius due to the effects described at the end of ^4.4.11 
The density profile of this test compares very well with that of the corresponding test of 
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Figure 4.13. Radial density profile of a collapse simulation when nested splitting was applied within 
a radius of 2 x 10^'^ pc, with particles initially taken from a settled distribution. The green points 
show fine particles and the black points coarse particles. The red line indicates the p cx profile. 
This is the end of the simulation at t ~ 1.01 tff. 

nesting splitting (i.e. nested splitting with a sink particle). This test extends the conclusion 
that the boundary effects due to on-the-fly splitting are not significant in a case where most 
coarse particles (~ 65%) are split. 

4.4.1.5. Isothermal collapse using nested splitting with particles initially taken 
from a settled distribution 

We have repeated the collapse simulation applying nested splitting to the centre of a cloud 
of 10,482 particles initially taken from a settled distribution with a uniform density of p = 
3.74 X 10^^^ g cm~'^. We have used the same initial conditions. We have simulated collapse 
with the same code as above f ^4. 4. 11 1). The fine region has the same radius of 2 x 10~^ pc. 
There is no sink particle. Fig. 14.131 shows the radial density profile of the cloud at t = 1.01 
tff. The density complies well to the predicted profile, p oc r^^. It compares very well with 
the simulation of ^4.4.11 1. when particles were initially taken on a lattice (top panel of Fig. 
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Figure 4.14. Initial state of a stable isothermal sphere. Top: Radial density profile of the isothermal 
sphere. The red line indicates the solution of Eqn. 12.541 Bottom : Thin equatorial slice (Az = 4 x 10^'^ 
pc) of the isothermal sphere showing the boundary between the coarse (black points) and the fine 
(green points) particles, immediately after nested splitting was applied within a radius of 3 x 10^^ pc. 
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14.111) . Note that in Fig. I4.13l there is more noise that in the top panel of Fig. 14.111 The cloud 
has contracted to a radius of 6.5 x 10~^ pc. The maximum density at the centre of the cloud 
is Ppeak = 1-8 X lO"^'^ g cm~'^. There are 74,394 particles in total, with 5,326 coarse particles 
having been split. This test rules out the possibility that the results presented above are 
biased due to some preferred orientation of the initial lattice. 

4.4.2 Stable Isothermal Sphere 

We test particle splitting on a simulation of a stable isothermal sphere. Both methods 
show again good results. The new calculation of h ( HA.'AA\i has clearly eliminated previous 
inefficiencies in treating the fine region boundary. Fine particles still move towards the gaps 
between coarse particles at the boundary, but there is little penetration through the boundary. 
In particular, as the fine region slowly expands, coarse particles trapped in the expanding fine 
region are automatically split. Any fiuctuations input by the application of particle splitting 
are being removed and the sphere settles towards the predicted density profile after a few 
free-fall times. 

The initial conditions consist of a spherical cloud of mass M = 1 Mq and radius R = 0.1 
pc. The cloud has central density Ppeak = 3.04 x lO"^'' g cm^'^ and uniform temperature T 
= 7.9 K. The density at the outer edge of the cloud is Pedge = 1-05 x 10^^'^ g cm^'^. The 
cloud is stable as it has H = 3 ( ^2.11.3|) . There is no rotation. We have used 10,185 particles 
to simulate the cloud, initially taken on a lattice. Their positions are calculated according 
to the scheme given in ^2.11.31 For this test we use our complete self-gravitating SPH code 
(chapter 121) . The top panel of Fig. 14.141 shows the initial radial density profile of the cloud. 
We use crosses as the particles initially are on a lattice and the radial profile therefore consists 
of very few different points. The red line indicates the solution of Eqn. 12.541 The bottom 
panel of Fig. 14.141 is a thin equatorial slice (Az = 4 x 10^'^ pc) of the isothermal sphere 
showing the boundary between the coarse (black points) and the fine (green points) particles, 
immediately after nested splitting was applied within a radius of 3 x 10~^ pc. 

We present three different simulations for the evolution of a stable isothermal sphere: 
using nested splitting, using on-the-fly splitting and using nested splitting with particles 
initially taken from a settled distribution. 
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4.4.2.1. Stable isothermal sphere using nested splitting 

The top panel of Fig. 14.151 shows the radial density profile of the isothermal sphere after 
t = 8.17 tff, when nested splitting was applied with the fine region having a radius of 3 x 
10-^ pc. The green points show fine particles and the black points coarse particles. The 
red line indicates the solution of Eqn. 12.541 The density profile compares well with the 
initial one (top panel of Fig. I4.14() . There are 22,833 particles in total with 1,054 coarse 
particles having been split. There are still some particles on either side of the boundary 
between the coarse and the fine regions having their density incorrectly modelled. These 
are basically new fine particles (split a few time-steps before) that have taken inappropriate 
initial values for their smoothing length. Their smoothing length and density will very soon 
settle to the expected values. However, the number of these particles is significantly reduced 
due to the application of the new method for calculating h. Moreover, the deviation from the 
expected density profile for these particles is greatly reduced (c/. with the top panel of Fig. 
14. 9j) . The introduction of new fine particles into the simulation creates the opposite effect to 
neighbouring coarse particles: their smoothing length is over-estimated for a few time-steps 
and thus their density is under-estimated. 

The fine particles inside the fine region boundary are still attracted by the coarse particles 
just outside the fine region. The fine region appears to expand slowly. However, there is now 
a clear boundary between the coarse and fine regions (bottom panel of Fig. 14.151 cf. with the 
bottom panel of Fig. 14. 9|) . Fine particles do penetrate through the coarse particles. However, 
the majority of the latter are always kept outside the fine region. In fact, as fine particles 
slowly move to a region of lower density some coarse particles move inside the fine region 
radius and are automatically split. The coarse particles surrounded by fine particles in the 
bottom panel of Fig. 14.151 are split in the subsequent time-steps. This way, the number of 
coarse particles being split increases from 459 initially to 1,054 at the end of the simulation. 
The voids in the fine region shown in the equatorial slice of the bottom panel of Fig. 14. 151 are 
projection effects produced when coarse particles exit the thin equatorial slice. Note that the 
thickness of the equatorial slice is less than, or of the order of the smoothing length of the 
coarse particles {Az = 4 x 10~'^ pc). Coarse particles exit temporarily the equatorial slice 
due to numerical noise in their inward motion (this motion is caused by the slowly expanding 
fine region). There is no way one can damp this noise at a scale smaller than the smoothing 
length of the particles. 
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Figure 4.15. Evolution of a stable isothermal sphere after t = 8.17 t//, when nested splitting was 
applied within a radius of 3 x 10^^ pc. The green points show fine particles and the black points coarse 
particles. Top: Radial density profile of the isothermal sphere. The red line indicates the solution of 
Eqn. 12.541 Bottom : Thin equatorial slice (Az — 4 x 10^^ pc) of the isothermal sphere showing the 
boundary between the coarse and the fine particles. 
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Figure 4.16. Evolution of a stable isothermal sphere after t ~ 8.59 tff, when on-the-fly splitting was 
applied. The density threshold is taken to be pmax = 2.48 x 10^'^" g cm"'^. The green points show fine 
particles and the black points coarse particles. Top: Radial density profile of the isothermal sphere. 
The red line indicates the solution of Eqn. 12.541 Bottom : Thin equatorial slice (Az = 4 x 10^"^ pc) 
of the isothermal sphere showing the boundary between the coarse and the fine particles. 
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Figure 4.17. Evolution of a stable isothermal sphere when nested splitting was applied within a 
radius of 3 x 10^^ pc, with particles initially taken from a settled distribution. The green points show 
fine particles and the black points coarse particles. The red line indicates the solution of Eqn. 12.541 
Radial density profile of the simulation after t — 6.63 t//. 

We conclude that SPH with nested splitting constrains the cloud to remain at the same 
overall equilibrium state at all times, while the fine region boundary is permitted to evolve 
slowly within this global profile. In fact, the particles whose density deviates form the ex- 
pected profile quickly have their density settled to the predicted values. This test shows that 
the boundary effects introduced by the application of nested splitting are not significant as 
the stable isothermal sphere retains its overall properties. It also demonstrates the efficiency 
of the new method for calculating h 

4.4.2.2. Stable isothermal sphere using on-the-fiy splitting 

The top panel of Fig. 14.161 shows the radial density profile of the isothermal sphere after t = 
8.59 tff, when on-the-fly splitting was applied with density threshold Pmax = 2.48 x 10"^*^ g 
cm~'^. The green points show fine particles and the black points coarse particles. The red line 
indicates the solution of Eqn. 12.541 The density profile compares very well with the initial 
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one (top panel of Fig. I4.14() . There are 27,057 particles in total with 1,406 coarse particles 
having been split. The density profile is similar to that of the nested splitting simulation 
(top panel of Fig. I4.15|) . Similar boundary effects can be observed in this test. There are 
some particles on either side of the boundary between the coarse and the fine regions having 
their density incorrectly modelled. These are basically new fine particles (split a few time- 
steps before) that have taken inappropriate initial values for their smoothing length. Their 
smoothing length and density will very soon settle to the expected values. The introduction 
of new fine particles into the simulation creates the opposite effect on neighbouring coarse 
particles: their smoothing length is over-estimated for a few time-steps and thus their density 
is under-estimated. 

The fine particles inside the fine region are attracted by the coarse particles just outside 
the fine region. The fine region appears to expand slowly. However, there is again a clear 
boundary between the coarse and fine regions (bottom panel of Fig. 14.161 cf. with the bottom 
panel of Fig. 14.91 note that in on-the-fly splitting the fine region's radius is larger). Fine 
particles do penetrate through the coarse particles. However, the majority of the latter are 
always kept outside the fine region. In fact, as fine particles slowly move to a region of 
lower density some coarse particles have their density estimates increased above the density 
threshold and are automatically split. The coarse particles surrounded by fine particles in 
the bottom panel of Fig. 14. 161 are split in the subsequent time-steps. This way, the number of 
coarse particles being split increases from 959 initially to 1,406 at the end of the simulation. 
Therefore, a smaller fraction of coarse particles is split compared to the nested splitting 
simulation. The voids in the fine region shown in the equatorial slice of the bottom panel of 
Fig. 14. 161 are again projection effects produced when coarse particles exit the thin equatorial 
slice. Note that the thickness of the equatorial slice is less than, or of the order of the 
smoothing length of the coarse particles {Az = 4 x 10^^ pc). As in the previous simulation, 
coarse particles exit temporarily the equatorial slice due to numerical noise in their inward 
motion. There is no way one can damp this noise at a scale smaller than the smoothing 
length of the particles. 

We conclude that SPH with on-the-fiy splitting constrains the cloud to remain at the 
same overall equilibrium state at all times, while the fine region boundary is permitted to 
evolve slowly within this global profile. In fact, the particles whose density deviates form 
the expected profile quickly have their density settled to the predicted values. This test 
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shows that the boundary effects introduced by the apphcation of on-the-fly sphtting are not 
significant as the stable isothermal sphere retains its overall properties. In this test, on-the-fly 
splitting produces similar results to nested splitting, again more efficiently. 

4.4.2.3. Stable isothermal sphere using nested sphtting with particles initially 
taken from a settled distribution 

We have repeated the simulation for the evolution of the stable isothermal sphere applying 
nested splitting to the centre of a cloud of 10,482 particles initially taken from a settled 
distribution of uniform density. The particles are then positioned according to the scheme 
given in H2.11.31 We have used the same initial conditions and the same code as above 
(' j4.4.2l l). The fine region has the same radius of 3 x 10~^ pc. Fig. 14.171 shows the radial 
density profile of the cloud at t = 6.63 tfj. The density complies well to the initial profile (top 
panel of Fig. I4.14|) . It compares very well with the simulation of H4.4.2I 1. when particles were 
initially taken on a lattice (top panel of Fig. I4.15j) . There are 22,602 particles in total, with 
1,010 coarse particles having been split. The boundary effects are similar to those discussed 
above. In the simulation of ^4.4.21 1. at t = 6.63 tjf, roughly the same number of coarse 
particles had been split. This test rules out the possibility that the results presented above 
are biased due to some preferred orientation of the initial lattice. 

4.4.3 Rotating Cloud 

Having tested the new method against a problem of homogeneous inflow as well as the 
evolution of a stable cloud, we now apply it to the standard test simulation, discussed in 
chapter 131 Application of the new method to a more realistic problem with known solution 
may reveal disadvantages of the new method that we have not calculated or thought of. 

We apply the initial conditions used in the simulations of the previous chapter f ij3.1() . 
i.e. we use uniform-density, isothermal (cq = 0.17 km s~^; T = 7.9 K; a ~ 0.26), rotating 
{n = 7.2 X 10-13 j-ad s-^ (3 ^ 0.16), spherical clouds of mass M = 1 Mq and radius R ~ 
0.02 pc, with particles cut initially from a face-centred cubic lattice and then given an m = 2 
azimuthal perturbation by adjusting their spherical polar azimuthal coordinate, (p, to a value 
(j)* given by 

<j) = ct>* + with A = 10% amplitude. 

m 

We use our complete self-gravitating SPH code with 50 neighbours (chapter [21), including 
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adiabatic heating. Our findings in the tests of the previous chapter, and in the tests of 
particle sphtting in the present chapter, assist in fine-tuning the values of the initial number of 
particles and the density above which adiabatic heating switches on. This way, the simulations 
always obey the Jeans condition with minimum computational cost. 

In chapter |3J we have shown that convergence to the results of Truelove et al. (|19971 
I1998|l and Klein et al. H1999|) can be achieved only with high resolution simulations. We have 
also shown that filamentary singularities can be obtained only if the simulations are evolved 
isothermally for sufficient time. The new method allows us to start a coarse simulation with 
low- to medium-resolution until it reaches the maximum resolvable density, Pmax (Eqn. 14.3)1 . 
when particle splitting can be applied to increase the resolution 13-fold. If adiabatic heating 
starts before the fine simulation reaches its resolution limit, then the simulation obeys the 
Jeans condition at all times. Specifically, in the adiabatic heating regime the Jeans mass 
increases with increasing density (c/. Eqns. IA.8I &: 0.39|) so that Eqn. 14.11 is always valid as 
long as adiabatic heating switches on before the fine simulation reaches its resolution limit. 
The higher the initial number of particles the longer the simulation can evolve isothermally, 
allowing for the use of a higher value for pQ, the density above which adiabatic heating starts. 

In order to use a high density for the switch to the adiabatic regime (po=10~^^ g cm~'^), 
we have to start the simulation with no less than 40,000 particles. For comparison with the 
simulation of ^3.3.1) we also present a simulation with pQ=10~^'^ g cm^'^. In order for such 
a simulation to obey the Jeans condition at all times, a cloud of as low as 10,000 particles 
initially can be used. In our effort to evolve both simulations with minimum computational 
cost, we use the minimum initial number of particles that allows the Jeans condition to be 
obeyed at all times for the above values of po- In particular, the simulation with pQ=10~^^ 
g cm~'^ starts with 40,000 particles and the simulation with /3o=10~^^ g cm^^ starts with 
10,000 particles. We apply both versions of the new method to each simulation. 

Bate & Burkert (|1997|1 have shown that clouds of 10,000 particles do not resolve the binary 
formation. In contrast, clouds of 40,000 particles do resolve the formation of the binary and, 
in general, produce results very similar to those obtained with isothermal simulations of 
higher resolution. 

The particle splitting simulation with 40,000 particles initially gives results consistent with 
those of Truelove et al. and Klein et al. as well as those presented in ^3.3.41 for the 600,000 
particle simulation, but with fewer particles, demonstrating the economy in computation 
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achieved with particle sphtting. The new method reduces the computational cost on two 
counts: we do not have to use high resolution from the beginning of the simulation and when 
we do increase the resolution, this does not have to happen everywhere but only in regions 
of interest. 

Using particle splitting just before the 10,000 particle simulation violates the Jeans con- 
dition, proves to be sufficient to reproduce results of higher resolution simulations (Bate &: 
Burkert (nW|) . gS3IIlE3I2l2), but not the results of Truelove et al. ipHTl [T!inH|l . Klein et 
al. (|1999)1 and ^8.3.41 This happens because the initial stages of the cloud evolution are not 
properly modelled due to low resolution. We present a simulation that cannot give results 
consistent with those of Truelove et al. and Klein et al. despite the fact that it obeys the 
Jeans condition at all times, to demonstrate that particle splitting is a necessary, but not a 
sufficient, condition for the reliability of a simulation. 

Both versions of the new method produce similar results for each simulation. On-the- 
fly splitting is again more efficient than nested splitting in terms of computational cost. 
Tables Wj\ fc 14.31 show a summary of our results. For comparison purposes, they also list the 
corresponding results of Truelove et al. H1998() and Klein et al. H1999|) . 

We do not need to present the simulation with particles initially taken from a settled 
distribution. The only case that the bar fragments is the nested splitting simulation for a 
cloud of 10,000 particles initially and the results of such a simulation with particles initially 
taken from a settled distribution are similar to those of ^3.61 Therefore, we refer to the 
discussion of the results of H3.6I 

All figures presented here are column density plots viewed along the rotation axis. The 
captions indicate the units of the colour tables. They also give the linear size of the figure 
and the time of the simulation. 

4.4.3.1. Nested Splitting for a Cloud of 40,000 Particles 

The top panel of Fig. 14.181 is a column density plot viewed along the rotation axis and shows 
the density projected on the x-y plane initially (at t = 0), where the density enhancements 
indicate the m=2 perturbation. In a medium to high resolution simulation of a rotating, 
spherical, uniform, isothermal cloud with an m = 2 perturbation, the binary separation is 
expected to be ~0.004 pc ( ^3.3|) . Taking into account the flattened shape that the cloud 
takes after ~1 tff, the fine region in nested splitting takes a cylindrical shape and its radius 
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is large enough to contain the binary formed and any other structure associated with it (e.g. 
spiral tails). The cylinder's height is as large as the thickness of the disc formed (~l/4 of the 
disc radius). Specifically, the fine region radius is 0.003 pc. 

Nested splitting is applied after t =1.244 tjf, when the binary has started forming and 
the maximum density is about to exceed the density threshold, Pmax = 2.4 x 10^"'^'^ g cm^'^, 
above which the simulation stops resolving the Jeans mass. Initially, 17,400 coarse particles 
lie inside the fine region radius and are split. The bottom panel of Fig. 14.181 is a column 
density plot immediately after the application of nested splitting. Note the smooth transition 
of the column density through the fine region boundary (its diameter is 0.006 pc while the 
linear size of the bottom panel of Fig. 14. 181 is 0.008 pc). 

In subsequent time-steps, all particles crossing the fine region radius are split on-the- 
fly. A well-defined binary forms and between the binary components a bar grows. The fine 
simulation would reach its resolution limit at pmax = 4 x 10"^^ g cm~^, but adiabatic heating 
starts at po = 10"^^ g cm~'^ and the simulation obeys the Jeans condition at all times. The 
top panel of Fig. 14.191 shows the column density just before adiabatic heating initiates at 
t =1.254 tff. 

The top panel of Fig. 14. 191 compares well with Fig. 12 of Truelove et al. (|1998|) : there are 
two elongated objects and a thin bar connecting them {ppeak = 7.14 x 10"^'^ g cm~^). 

The bottom panel of Fig. 14.191 shows the column density at the end of the simulation, 
t =1.265 tff. There are 281,331 particles in total with 20,020 coarse particles having been 
split. The mass of each fragment is 0.02 M© and its radius ~7 AU. Their separation is 515 
AU. The peak density of the simulation has risen to Ppeak = 1-97 x 10"^*^ g cm~^. The 
bar has not fragmented. There is a well-defined contrast between the density of the binary 
components and the bar that connects them (2 orders of magnitude). 

The final results compare well with those of Klein et al. (11999(1 (their Fig. 2 - note that 
in Klein et al. adiabatic heating starts about one order of magnitude earlier). 

The resolution of the simulation is sufficient that we can trust the result that the bar does 
not fragment up to t =1.265 tff. In fact, the resolution of this simulation is equivalent to the 
resolution of a 520,000 particle simulation and the results are indeed similar to those of the 
simulation with the highest resolution we conducted without splitting in ^3.3.41 f600.000). In 
particular, the peak density of the simulation is comparable to that of the simulation with 
600,000 particles (see discussion in ^3.3.4j) . This indicates the efficiency of nested splitting as 
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Figure 4.18. Nested Splitting for a Cloud of 40,000 Particles: Column density plots of the initial 
sphere and the cloud after splitting is applied. Top: Column density plot of the cloud initially {t = 
0). The linear size of this plot is 0.04 pc. The colour table has units of 1.18 x 10^ g cm~^. Bottom : 

Column density plot of the cloud after splitting is applied {t = 1.244 tff). The linear size of this plot 
is 0.008 pc. The colour table has units of 2.95 x 10'' g cm~^. 
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Figure 4.19. Nested Splitting for a Cloud of 40,000 Particles: Column density plots of the cloud 
before heating is applied and at the end. Top: Column density plot of the cloud before heating is 
applied {t = 1.254 tff). The linear size of this plot is 0.004 pc. The colour table has units of 1.18 x 

10® g cm^^. Bottom : Cohimn density plot of the cloud at the end (t — 1.265 tff). The linear size of 
this plot is 0.004 pc. The colour table has units of 1.18 x 10® g cm~^. 
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Figure 4.20. On-the-fly Splitting for a Cloud of 40,000 Particles: Column density plots after splitting 
is applied and before heating starts. Top: Column density plot of the cloud after splitting is applied 
{t = 1.251 tff). The linear size of this plot is 0.008 pc. The colour table has units of 2.95 x 10^ g 
cm~^. Bottom : Column density plot of the cloud before heating starts (t = 1.258 tff). The linear 
size of this plot is 0.004 pc. The colour table has units of 1.18 x 10* g cm~^. 
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Figure 4.21. On-thc-fly Splitting for a Cloud of 40,000 Particles: Column density plot of the cloud 
at the end {t — 1.277 i//)- The linear size of this plot is 0.004 pc. The colour table has units of 1.18 
X 10^ g cm~^. 

only 281,331 particles have been used. 

4.4.3.2. On-the-fly Splitting for a Cloud of 40,000 Particles 

On-the-fly splitting is applied after t =1.244 tff, when the maximum density is about to 
exceed the density threshold, pmax = 2.4 x 10~^^ g cm~^, over which the simulation stops 
resolving the Jeans mass. Initially, 2,955 coarse particles are split. The top panel of Fig. 14.201 
is the first column density plot after the application of on-the-fly splitting (t= 1.251 t//)- 
Note the smooth transition of the column density through the fine region boundary (within 
the red coloured area). 

In subsequent time-steps, all particles whose density exceeds the density threshold are split 
on-the-fly. A well-defined binary forms and between the binary components a bar grows. The 
fine simulation would reach its resolution limit at Pmax = 4 x 10~^^ g cm~^, but adiabatic 
heating starts at po = 10~^^ g cm~^ and the simulation obeys the Jeans condition at all 
times. The bottom panel of Fig. l4.2Ul shows the column density just before adiabatic heating 
starts at t =1.258 tff. 

The bottom panel of Fig. 14. 2UI compares well with Fig. 12 of Truelove et al. (|1998|1 : there 
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are two elongated objects and a thin bar connecting them {ppeak = 6.86 x 10"^'^ g cm~^). 

Fig. 14.211 shows the column density at the end of the simulation, t =1.277 tfj. There 
are 150,135 particles in total with 9,087 coarse particles having been split. The mass of each 
fragment is 0.03 Mq and its radius ~10 AU. Their separation is 475 AU. The peak density 
of the simulation has risen to Ppeak = 1-85 x lO^^'^ g cm~^. The bar has not fragmented. 
There is a well-defined contrast between the density of the binary components and the bar 
that connects them (2 orders of magnitude). 

The final results compare well with those of Klein et al. ()1999|) (their Fig. 2 - note that 
in Klein et al. adiabatic heating starts about one order of magnitude earlier). 

The resolution of the simulation is sufficient that we can trust the result that the bar does 
not fragment up to t =1.277 tjf. In fact, the resolution of this simulation is equivalent to the 
resolution of a 520,000 particle simulation and the results are indeed similar to those of the 
simulation with the highest resolution we conducted without splitting in ^233(600,000). In 
particular, the peak density of the simulation is comparable to that of the simulation with 
600,000 particles (see discussion in H3.3.4() . This indicates the efficiency of on-the-fly splitting 
as only 150,135 particles have been used. 

Comparison between the results of the on-the-fly splitting and the nested splitting sim- 
ulations clearly indicates that the results are similar (see Tables 14.21 & 14. 3() , although the 
on-the-fly splitting simulation has progressed more in time (the binary components have 
accreted more matter, they are closer and the bar is more inclined). However, on-the-fly 
splitting is much more efficient in obtaining them (using 131,196 less particles in total). 

4.4.3.3. Nested Splitting for a Cloud of 10,000 Particles 

The top panel of Fig. 14.221 is a column density plot viewed along the rotation axis and shows 
the density projected on the x-y plane initially (at t = 0) for a cloud of 10,000 particles, 
where the density enhancements indicate the m=2 perturbation. The coarser resolution used 
here does not change the initial appearance of the cloud (c/. the top panel of Fig. 14.18(1 . 
Following the same reasoning as in ^4.4.31 1. the fine region is cylindrical with the radius now 
taken to be 0.002 pc (the expected bar length is extracted from Fig. 7 of Bate & Burkert 

(HMD). 

Nested splitting is applied after t =1.105 t//, when the maximum density is about to 
exceed the density threshold, pmax = 1-5 x 10^^^ g cm~^, over which the simulation stops 
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Figure 4.22. Nested Splitting for a Cloud of 10,000 Particles: Column density plots of the initial 
sphere and the cloud after splitting is applied. Top: Column density plot of the cloud initially {t = 
0). The linear size of this plot is 0.04 pc. The colour table has units of 1.18 x 10® g cm~^. Bottom : 

Column density plot of the eloud after splitting is applied {t = 1.105 <//). The linear size of this plot 
is 0.008 pc. The colour table has units of 2.95 x 10'' g cm~^. 
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Figure 4.23. Nested Splitting for a Cloud of 10,000 Particles: Column density plots of the cloud 
before heating is applied and at the end. Top: Column density plot of the cloud before heating is 
applied {t = 1.265 t/j)- The linear size of this plot is 0.004 pc. The colour table has units of 1.18 x 

10^ g cm~^. Bottom : Column density plot of the cloud at the end (t = 1.301 t//). The linear size of 
this plot is 0.004 pc. The colour table has units of 1.18 x 10* g cm~^. 
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resolving the Jeans mass. At this time, the binary has not formed yet and the centre of 
the cloud has just started collapsing after the initial expansion phase ({Bate &: Burkert 1997|1 . 
Initially, 705 coarse particles lie inside the fine region radius and are split. The bottom panel 
of Fig. 14.221 is a column density plot immediately after the application of nested splitting. 
Note that the lattice which the particles were initially taken on has not broken yet. Also note 
the smooth transition of the column density through the fine region boundary (its diameter 
is 0.004 pc while the size of the bottom panel of Fig. 14.221 is 0.008 pc). The density of the 
fine particles is emphasised due to effects similar to those discussed in ^4.4.21 1 and the small 
dynamic range of the bottom panel of Fig. 14.221 

In subsequent time-steps, all particles crossing the fine region radius are split on-the-fiy. 
A central density enhancement grows and later forms a bar. On either side of the bar spiral 
tails form due to rotation. At both ends of the bar two objects finally form at t = 1.265 
tff. The fine simulation would reach its resolution limit at Pmax = 2.53 x 10"^'^ g cm~^, but 
adiabatic heating starts at po = 10^^^ g cm~^ and the simulation obeys the Jeans condition 
at all times. The top panel of Fig. 14.231 shows the column density just before adiabatic 
heating starts at t =1.265 tff. Subsequently, the bar fragments producing three fragments, 
one at the centre and two at equal distances from the centre. 

This is in accordance with the findings of the low resolution run of Bate & Burkert p997jl , 
where the binary formation is not resolved and just a bar forms that later fragments into 
multiple fragments. However, in this simulation, with application of nested splitting, the 
binary forms, although at a later stage than expected. 

The bottom panel of Fig. 14.231 shows the column density at the end of the simulation, 
t =1.301 tff. There are 55,125 particles in total with 3,794 coarse particles having been split. 
The mass of the binary fragments is 0.06 M0 and of the bar fragments 0.01 M0. The binary 
fragments have radii of ~25 AU. The bar is 410 AU long. The peak density of the simulation 
has risen to Ppeak = 2.14 x 10^^^ g cm~'^. The bar fragments as it does in the corresponding 
runs of Bate & Burkert (fTWTfl and those of ^^TU l3X2l 2. 

Clearly, the evolution of this simulation does not imitate the evolution predicted by the 
high resolution simulations of Truelove et al. (jTWlIT^^ . Klein et al. ((TO^ and ^3.3.41 (see 
Tables 10 
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Figure 4.24. Simulation of a Cloud of 130,000 Particles without Particle Splitting: Column density 
plot of the cloud at the end {t =1.291 tff). The linear size of this plot is 0.004 pc. The colour table 
has units of 1.18 x 10^ g cm~^. 

The equivalent resolution from the beginning of the simulation At the centre of this 
simulation the resolution is equivalent to that of a simulation with 130,000 coarse particles. 
We have constructed such a test to compare with the above simulation. We have used the 
same initial conditions with 130,000 particles initially taken on a lattice. Adiabatic heating 
initiates at po = 10^ ^'^ g cm~^. Particle splitting is not applied. 

Fig. 14.241 shows the end of the 130,000 particles simulation, at t =1.291 tff. In this 
case, the binary formation is resolved (forms at t =1.239 tff). The bar between the binary 
components does fragment, i.e. the results are similar to the simulations of Bate & Burkert 
(Pn7|l and those of EH I33I2I2. 

The difference between the results of this run and the results of the above simulation with 
only 10,000 particles initially is due to the low resolution initial phase of the latter simulation. 
In particular, in the latter simulation the initial expansion and the subsequent collapse that 
leads to the formation of the binary ([Bate fc Burkert 1997|) are not properly modelled (due 
to poor sampling) and this causes the delayed formation of the binary. The application of 
nested splitting and the fact that the Jeans condition is obeyed clearly assists in the eventual 
formation of the binary (c/. the result of the 10,000 particle simulation of Bate & Burkert 
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Figure 4.25. On-the-fly Splitting for a Cloud of 10,000 Particles: Column density plots of the cloud 
after splitting is applied and before heating starts. Top: Column density plot of the cloud after 
splitting is applied {t = 1.142 t//). The linear size of this plot is 0.008 pc. The colour table has units 
of 2.95 X 10^ g cm^^. Bottom : Column density plot of the cloud before heating starts {t ~ 1.293 i//). 
The linear size of this plot is 0.004 pc. The colour table has units of 1.18 x 10® g cm~^. 



4.4. TESTS 



121 



The figure is provided separately 
due to its large size 



Figure 4.26. On-tlie-fly Splitting for a Cloud of 10,000 Particles: Column density plot of the cloud 
at the end {t = 1.326 i//)- The linear size of this plot is 0.004 pc. The colour table has units of 1.18 
X 10^ g cm"^. 

(Unnii))- 

This conclusion is confirmed by the results of simulations of 20,000 and 30,000 particles 
initially, with application of on-the-fly splitting. In particular, we find that in these two runs 
the binary formation is resolved. Moreover, the evolution of the 30,000 particle simulation is 
closer to the converged solution of Bate & Burkert (|1997j) . However, application of particle 
splitting, despite assisting this simulation in creating the binary that would not be formed 
otherwise, can not guarantee solving all problems, due to poor sampling. 

4.4.3.4. On-the-fly Splitting for a Cloud of 10,000 Particles 

On-the-fly splitting is applied to the cloud of 10,000 particles after t =1.105 tff, when the 
maximum density is about to exceed the density threshold, Pmax = 1-5 x 10"^^ g cm~^, over 
which the simulation stops resolving the Jeans mass. Initially, 3,146 coarse particles are split. 
The top panel of Fig. 14.251 is the first column density plot after the application of on-the-fly 
splitting (t=1.142 tff). Note the smooth transition of the column density through the fine 
region boundary (within the green coloured area). Note that the lattice which the particles 
were initially taken on has not broken yet. Within the few time-steps elapsed after the first 
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coarse particles were split, a rugged appearance similar to the bottom panel of Fig. 14.221 
has been smoothed out as it did in ^4.4.21 2. Due to the low density threshold, many coarse 
particles are split on the fly within these few time-steps. At t=1.142 tjf, the binary has not 
formed yet and the centre of the cloud has just started collapsing after the initial expansion 
phase (|Bate Burkert 1997|1 . 

In subsequent time-steps, all particles whose density exceeds the density threshold are 
split on-the-fly. A central density enhancement grows and later forms a bar. On either side 
of the bar spiral tails form due to rotation. The fine simulation would reach its resolution 
limit at Pmax = 2.53 x 10~^'^ g cm~^, but adiabatic heating starts at po = 10"^'^ g cm~^ and 
the simulation obeys the Jeans condition at all times. 

The bottom panel of Fig. 14.251 shows the column density just before adiabatic heating 
starts at t =1.293 tjf. At both ends of the bar two objects finally form at t = 1.317 tjf. 
Subsequently, the bar does not fragment. This result is in accordance with the findings of 
the high resolution simulations of Truelove et al. (|TW1 [T!inH|l . Klein et al. I^HM^ and ^3.3.41 
where the bar does not fragment. However, in this simulation, due to the initial coarse 
resolution phase, the binary forms at a much later stage than expected. Fig. 14.261 shows the 
column density at the end of the simulation (t =1.326 tfj). The density contrast between the 
binary components and the bar that connects them is greater than one order of magnitude. 
There are 96,255 particles in total with 7,219 coarse particles having been split. The mass of 
each binary fragment is 0.06 Mq and its radius ~40 AU. The bar is 270 AU long. The peak 
density of the simulation has risen to Ppeak = 5.04 x 10~^^ g cm~'^. 

In this simulation the initial expansion and the subsequent collapse that leads to the 
formation of the binary ()Bate &: Burkert 1997|1 are not properly modelled (due to poor sam- 
pling) and this causes the delayed formation of the binary. The application of on-the-fly 
splitting and the fact the Jeans condition is obeyed clearly assists in the eventual formation 
of the binary (c/. the result of the 10,000 particle simulation of Bate & Burkert p997)l ). The 
fact that the binary forms later than in all other low resolution simulations of ij3.4l as well 
as the above simulation (nested splitting with 10,000 particles initially) can be attributed to 
the large number of coarse particles being split almost simultaneously at the initial stages. 
During the collapse phase following the initial expansion, many coarse particles obtain den- 
sity larger than the low density threshold we have used. We have used such a low value for 
Pmax (1-5 X 10^^^ g cm^'^) due to the small initial number of particles. The simultaneous 
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splitting of few thousand particles introduced excess density perturbations to the simulation. 
SPH and the new method for calculating h have been shown to be efficient in eliminating 
such perturbations. This is why we finally obtain the right result. 

Comparison between the results of the on-the-fly splitting and the nested splitting simula- 
tions with 10,000 particles initially, shows that only with on-the-fly splitting we can reproduce 
the results of the high resolution simulations of Truelove et al. (|1997| I1998|l , Klein et al. (|1999|1 
and N3.3.4I (see Tables &: I4.3)l . However, with the low density threshold set due to the 
small initial number of particles the computational cost of the on-the-fly splitting simulation 
grows considerably. Therefore, accuracy has been achieved at the expense of computational 
efficiency. This is against one of the primary aims of the new method that is designed to 
gain both in accuracy and computational efficiency. This confirms our previous conclusion 
that particle splitting, the way it has been formulated in this work, is a necessary, but not a 
sufficient, condition for the reliability and efficiency of a simulation. 

4.4.4 Conclusions 

The first two tests f ^4.4.1l fc l4.4.2() have proven that the new method for calculating h f ^4.3.4|) 
has greatly reduced the boundary perturbations introduced by the application of particle 
splitting to the regions where fine and coarse particles are in contact. In particular, the 
collapse simulations f ^4.4.1() show that there is no outward propagation of the fine region 
boundary in cases of infiow, which are the flows usually investigated in simulations of Star 
Formation. It also shows that the fine region reproduces the expected density profiles. 

The simulation of the evolution of a stable isothermal sphere ( H4.4.2|) shows that there is 
a distinct boundary between the fine and the coarse regions. The boundary may expand in 
cases of quiescent evolution of stable spheres, but the overall density profile accurately obeys 
the predicted one, irrespectively of the size of the fine region. Both tests have revealed that 
the on the fly version of particle splitting is much more efficient in terms of computational 
cost. In principle, nested splitting could be superior to on-the-fiy splitting if the boundary 
perturbations mentioned in ^4.3.41 persisted after the implementation of the new method for 
calculating h, as in nested splitting the fine region is larger than it really needs to be. This 
could prevent the boundary effects from propagating inwards and corrupting the simulation. 
However, both tests show that the boundary effects of both versions of the new method 
have been eliminated to a large extent and that the insignificant errors produced by particle 
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Table 4.2. Summary of results for the standard test simulations using particle splitting (first part). For each simulation the second column gives 
the initial number of particles and the third column lists the density above which adiabatic heating starts. Pmax is the density at which the coarse 
simulation reaches its resolution limit. j, is the peak density at the end of the isothermal regime, tspi is the time at which particle splitting is 
applied, theat is the time at which heating starts and tt,i„ is the time the binary forms. The last two rows list the relevant results of the finite difference 
simulations of Truelove et al. I)1998|l and Klein et al. fll999|l . For these two rows, we have used the initial number of cells instead of the number of 
particles in column 2. 
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Table 4.3. Summary of results for the standard test simulations using particle splitting (second part). For each simulation the second column gives 
the initial number of particles and the last column the final number of particles. The third column lists the density above which adiabatic heating 
starts. is the peak density at the end of the simulation, t^nd is the time at the end of the simulation and ifrag is the time the bar fragments. 

The fifth column gives the number of fragments produced in the bar. The last two rows list the relevant results of the finite difference simulations of 
Truelove et al. H1998|l and Klein et al. H1999fl . For these two rows, we have used the initial and the final number of cells instead of the number of 
particles in columns 2 and 13 respectively. 
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splitting are similar in both versions of the method. Therefore, on-the-fly splitting can be 
safely used as it is superior in terms of computational cost. 

We have also applied particle splitting to collapse simulations like those of the previous 
chapter. We conclude that the simulations with 40,000 particles initially, that use a suffi- 
ciently high density for the switch to adiabatic heating {pQ=W^^^ g cm~'^) while obeying the 
Jeans condition at all times, reproduce the results of Truelove et al. (|1998|) . of Klein et al. 
(fTTM|l . and of our highest resolution simulation of gS321 (600,000), but with only ~200,000 
particles in total at the final stages. The resolution of these simulations is equivalent to the 
resolution of a simulation with 520,000 particles, which clearly indicates the efficiency of par- 
ticle splitting. In particular, with particle splitting we achieved ~ 60% economy in memory 
and ~ 70 — 75% economy in CPU. A summary of our results is given in Tables fc 14.31 

With particle splitting the simulation obeys the Jeans condition at all times with great 
computational gain. Our new method has succeeded in meeting the objectives we have set 
in i j4.21 We can now apply it to simulations of clump-clump collisions. 

On-the-fly splitting is more efficient than nested splitting, as even less particles are re- 
quired. This is due to the fact that in on-the-fly splitting no particles are split unnecessarily. 
On-the-fly splitting is our preferred version of the new method. We will use on-the-fly splitting 
in simulations of clump-clump collisions in the next chapter. 

In all simulations presented here, the Jeans condition is always obeyed and therefore there 
are no artificial fragments. This demonstrates the Jeans condition for fragmentation provides 
a very strong test for the significance of numerical results. 

However, with the simulations with 10,000 initially and po=10~^'^ g cm~^, we have shown 
that particle splitting in response to the imminent violation of the Jeans condition is a 
necessary, but not a sufficient, condition for the reliability and efficiency of a simulation. In 
particular, the low resolution of the initial stages of the cloud evolution has prevented the 
nested splitting simulation from reproducing the result of finite difference simulations of the 
same problem. It has also prevented the on-the-fly splitting simulation from obtaining this 
result efficiently. 



Chapter 5 



On-the-fly Splitting Simulations of 
Clump-Clump Collisions 



Simulations of cloud-cloud collisions ( Chapman et al. 1992} Pongracic et al. 19921 iTiirner et al. 19951 



IWhitworth et al. 1995| IBhattal et al. 1998|l evolve over several orders of magnitude in den- 
sity and involve fragmentation to produce protostellar objects. Therefore, it is essential that 
fragmentation is properly modelled, i.e. the resolution is sufficient that the Jeans condition 
is always obeyed. The particle splitting simulations of the standard test with only 10,000 
particles initially f ^4.4.3() have shown that the Jeans condition is a necessary but not suffi- 
cient condition for reliable simulations; in that case the initial expansion and the subsequent 
collapse of the cloud were not properly modelled. 

The cloud-cloud collision simulations of Bhattal et al. (|1998|1 were conducted with small 
numbers of particles (2,000-8,000 per cloud) and implemented with a constant gravity soft- 
ening, e. However, small numbers of particles lead to violation of the Jeans condition and 
inaccurate modelhng of the density field (ESSl El & SHJ)- Also Bate & Burkert ((TM7|l 
have demonstrated that the implementation of constant e-softening can be very misleading 
for SPH calculations with gravity, as it can both artificially induce fragmentation of Jeans 
stable lumps and/or stabilise Jeans unstable lumps against collapse, depending on the e/h 
ratio (SIJ- 

Therefore, it is worth conducting new high resolution simulations in which the Jeans 
condition is explicitly monitored and obeyed, in order to obtain a more detailed understanding 
of the processes involved. In this chapter, we perform such high resolution simulations of 
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clump-clump collisions by applying particle splitting, and in particular the on-the-fly splitting 
version. 

In chapter [IJ we have shown that on-the-fly splitting is both a reliable and a computa- 
tionally efficient method. Therefore, with on-the-fly splitting, we can overcome the problem 
of insufficient resolution for modelling fragmentation and so we can obtain credible results. 
We use /i-softening, i.e. we set e = h. All particles have an individual time-varying h. 

We perform three simulations starting from the initial conditions used by Bhattal et al. 
()1998|) in order to make direct comparisons with their results, and hence to quantify the effects 
of low resolution and e-softening on those results ( ^5.2j) . We also investigate fragmentation 
induced by collisions between low-mass clumps by performing a series of simulations with 
various combinations of relative velocity and impact parameter f ^5.3|) . We note that collisions 
of clumps from the low-mass end of the clump mass function have not been investigated before. 

Previous collision simulations have shown that all fragmentation events happen in the 
shocked gas layer formed at the interface of the collision. It is therefore very important 
that the formation of this layer is properly modelled. In the following section, we define 
the appropriate minimum initial number of particles for our simulations to obey the Jeans 
condition at all times, and to resolve the formation of this shock. We also define the physical 
and numerical initial conditions for our simulations. 



5.1 Initial conditions 

When two clouds collide, a layer of shocked gas forms at the interface between the clouds 

( Chapman et al. 1992||Pongracic et al. 1992tlT;7rner e.t al. 1995UWhitworth et al. 1995llRhatta,] e.t al 1998jl . 



The layer eventually fragments producing groups of protostars. Whitworth et al. ()1994b|) 
have shown that the mass of a fragment, Mjrag, produced in such a shocked layer is given by 



a? 



where is the effective sound speed of the gas in the shocked layer, pdoud is the initial 
pre-shock density of each one of the two clouds involved in the collision and M is the Mach 
number of the collision, defined as the ratio of the relative velocity of the clouds, VcoU, to 
{vcoii = "^v cloud-, where Vdoud is the bulk velocity of each cloud). 
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Table 5.1. Minimum initial numbers of particles ^)vos ^ -^pl^ " different values of Mq and M. so that the Jeans condition is obeyed at all times. 
In all cases, we have used -/V„=50. The fifth column lists the mass of the fragments produced in the layer. 
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The pre-shock density of each cloud, Pdoud, is given by 

Pdoud ~ ^3^2 ' (5-2) 

if we make the assumption that the clouds are in hydrostatic balance (i.e. in an equilibrium 
state). Mq is the cloud mass and ao is the effective pre-shock sound speed, i.e. a sound speed 
which represents thermal and turbulent pressure. 
Combining Eqns. IS.ll fc l^T^ we obtain 



(5.3) 



In order for the formation of these fragments to be resolved, their mass must be larger, or 
of the order of, the minimum mass resolvable by SPH, 2Mmin- Substituting Mmin = Nnmptci 
we obtain, 



Mq 
Ntotal 



(5.4) 



and hence 



(5.5) 



Using the scaling relation, 



ao ~ 0.2 km s"M — ^ 
(c/. Larson (|1981j) ). we finally obtain 

Ntotai>2N^(^^y M'/', (5.6) 

where we have used 0^=0.2 km s^^ for the post-shock sound speed, or equivalently rs=10 K 
for the post-shock temperature. This is a reasonable assumption since for typical post-shock 
densities, the gas cools very rapidly and efficiently down to 10 K. The third column of Table 
15. II lists the minimum initial numbers of particles, N^^^g, needed to resolve fragmentation in 
the shocked layer for the different values of Mq and M. In ^2.11.11 we have shown that the 
shock is adequately modelled with our code for Vcioud=^-^ km s~^ (7W=5). 

Above pq = 10"^^ g cm^^ the gas heats up adiabatically, as described in chapter [21 In 
chapter |31 we have used different values for pQ in order to explore the performance of our 
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code under different circumstances. Tlie value of po above which adiabatic heating switches 
on is defined by where heating due to gravitational collapse exceeds the rate of cooling by 
dust - because the dust becomes optically thick to its own radiation (e.g. Whitworth et al. 
JMHl))- po is given by 

m3/2(fcT)V2 
PO J^2 ' 

where h and k are the Planck and Boltzmann constants respectively and c is the speed of light. 
If we substitute the values of the constants, T=10 K and assume that the gas is molecular, 
then we arrive to a value for po ~ 10^^^ g cm~^. 

In order for the Jeans condition to be obeyed not just for the first fragments formed in 
the layer but also up to po = 10^^^ g cm~^, we need 

12iV„MoG3/2pV2 

^^otal > . (5.7) 

(Eqn. 14.2(1 . To find the minimum initial number of particles, A^™", needed for the Jeans 
condition to be obeyed in the layer up to po = 10~^^ g cm~^, we divide by 13 the number 
of particles given by Eqn. 14. 2[ as we can increase the number of particles 13-fold by using 
particle splitting once. The fourth column of Table 15.11 lists the values of -/V™™ for the 
different values of Mq. 

From Table IKTTl we observe that N'^^^ is always larger than N^^^g so that we always use 
-^po™ ^s the initial number of particles in our simulations. In all cases, we use A'^„=50. 

Collisions between low-mass clumps (Mo=10 Mq) seem to be the most computationally 
efficient to explore as they only need ~15,000 particles initially for the Jeans condition to 
be obeyed at all times. Collisions between such clumps have not been investigated before. 
In all cases, two equal mass clumps of 10 collide with opposite bulk velocity vectors. In 
molecular clouds with hierarchical substructure, collisions between clumps at the same level 
of the hierarchy (i.e. with comparable mass) are the most probable ()Scalo 1985|1 . 

We have conducted simulations with A^=5, 10, and 15 and with 6, the ratio of the impact 
parameter to the radius of the cloud, satisfying 6=0.0, 0.2, 0.4, 0.6 and 0.8. Whitworth et al. 
(|1994aj) have argued that the equations used to estimate Mfrag are approximately valid not 
only for head-on collisions (6=0.0) but also for collisions with 6 <0.7. The 6=0.8 collisions are 
conducted only for completeness, as the clumps do not interact strongly, even in the slowest 
collision (A^=5). 



132 



CHAPTER 5. CLUMP-CLUMP COLLISIONS 



We also repeat three simulations using the initial conditions of Bhattal et al. H1998|) with 
the Jeans condition obeyed at all times (seventh row in Table 15. If) . Our aim is to make 
a direct comparison of our results with theirs, in order to be able to quantify the effect of 
fragmentation being properly modelled and the benefit of using particle splitting. We perform 
the simulations with 6=0.2, 0.4 and 0.5, as the first and the third ones are categorised by 
Bhattal et al. as having different mechanisms in play for the formation of binary /multiple 
protostellar systems (rotational instability in a disc vs. fragmentation of the layer). We also 
perform the second one as Bhattal et al. claim that it combines both mechanisms. 

From Table ISTTl we see that for 10 Mq clumps, the pre-shock sound speed ao ~0.35 km 
s~^ so that the pre-shock temperature of the clumps is Tq ~35 K. For typical post-shock 
densities cooling of the gas is very rapid and efficient so that it cools down to 10 K in the 
shocked layer very quickly. In order to be able to model this behaviour, we have slightly 
altered our barotropic equation of state (Eqn. 12. 4j) so that 



P(r) 



1+ (pW 
' po 



4/3- 



p ^ Pi; 

1/2 (5.8) 
, P>Pi, 



where oq ~0.35 km s^^, ~0.2 km s~^ and pi=1.2 x 10^^'^ g cm^^. pi is the density 
at which cooling initiates (jTohline 1982|1 . This equation of state can be divided into four 
regimes. An isothermal regime with To=35 K for densities below pi. A cooling stage, ap- 
proximating to T oc /9~^/^ for densities just above pi, until the temperature reaches Ts=lQ 
K. The temperature then remains constant at Ts=10 K for densities below pQ. Finally, the 
gas heats up adiabatically above pQ = 10~^^ g cm~^. 

We can summarise our initial conditions as follows. Two equal mass clumps collide with 
equal but opposite velocities. In repeating the simulations of Bhattal et al. we use Mo=75 
Mq, oq ~0.6 km s ^ (corresponding to ro=100 K); A4=9 (corresponding to Vdoud— 

0.9 km 

s~^), 6=0.2, 0.4 and 0.5, and 110,000 particles per clump. In the second set of simulations, 
which we call low-mass clump collisions, we use Mo=10 Mq, ~0.35 km s""*^ (corresponding 
to To=35 K); M.=5, 10, 15 (corresponding to Vcioud=0.5, 1.0 and 1.5 km s^^), 6=0.0, 0.2, 
0.4, 0.6 and 0.8, and 15,000 particles per clump. 

In all cases, the clumps are stable isothermal spheres made with a procedure similar to 
that of ii2.11.3l The particles, before they are moved to their final positions, are relaxed 
to uniform density, as they were initially taken from a random distribution (this relaxed 
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The figure is provided separately 
due to its large size 



The figure is provided separately 
due to its large size 



Figure 5.1. Column density plots for the 6=0.2 collision of two 75 M0 clumps at the end of the 
simulation, t ^^0.64 Myr. Left : Column density plot viewed along the z-axis. The linear size of this 
plot is 0.04 pc. The colour table has units of 1.18 x 10^ g cm^^. Right: Column density plot viewed 
along the y-axis. The linear size of this plot is 0.028 pc. The colour table has units of 2.41 x 10^ g 

_2 

cm . 

distribution is what we called in previous chapters a "settled" distribution). The colliding 
clumps are taken to be touching before the simulations start, to save computational time. 
During the preceding approach of the colliding clouds, mutual tidal distortion will be small 
because 7W > 5. We use our full self-gravitating SPH code (chapter |2j with the above 
barotropic equation of state (Eqn. 15.8(1 : we use /i-softening with our TCG method ( H2.5\i : 
we use A'^„=50; and we use on-the-fly splitting at pmax = /5o/13^ = 6.0 x 10"^'' g cm~^, as 
at this stage the Jeans condition would otherwise stop being obeyed. After application of 
particle splitting, adiabatic heating starts before the fine simulation reaches its resolution 
limit. Therefore, the Jeans condition is obeyed all the way up to the highest densities we can 
achieve. 

We use false-colour column density plots to present our results. All structure formed is 
contained within a single layer and therefore such plots are not greatly confused by projec- 
tion effects. Column density plots are preferred to particle plots as the former give a more 
accurate representation of the total density field, and of what would be seen in optically thin 
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The figure is provided separately 
due to its large size 



Figure 5.2. Column density plot for the 6=0.2 collision of two 75 M© clumps viewed along the y-axis, 
at the end of the simulation, t -^0.64 Myr, showing the network of filaments produced. The linear size 
of this plot is 0.2 pc. The colour table has units of 4.73 x 10^ g cm~^. 



molecular-line (or dust-continuum) radiation, assuming a uniform excitation temperature (or 
dust temperature). We do not use contour plots, as we believe that colour column density 
plots are easier to read and give more information. The figure captions give the linear size of 
each plot and the plane onto which mass is projected. They also give the time of the snapshot 
they represent and the units of the colour table. We have taken the clouds always to collide 
along the x-axis, with offset parallel to the y-axis. 



In all cases, when a fragment forms, its linear density profile appears like a normal dis- 
tribution around the peak density. At the end of our simulations, this profile is very steep. 
To infer its mass, we take the fragment to extend out to 3-a, where a is the FWHM, and we 
measure the mass within this radius. This method gives good results for the mass and the 
radius of spherical as well as disc-like objects. 
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5.2 Repeating the Simulations of Bhattal et al. 
5.2.1 Simulation with Mo=75 M© and 6=0.2 

Bhattal et al. H1998|l found that for low-6 cohisfons a single spherical protostellar object 
formed at the centre of the collision and then accreted material from the spindle-shaped fila- 
ment. The protostellar object soon became disc-like and grew in mass and angular momen- 
tum as the offset between the opposing accretion flows along the spindle increased. Several 
accretion-induced rotational instability events took place as the primary protostellar disc 
ejected material and angular momentum through spiral arms. The ejected material inter- 
acted with the accretion flows and produced multiple lower mass companions to the massive 
primary. The primary had mass between 20-60 M©. The low-mass companions formed were 
not always resolved properly, as in some cases they contained less than 50 particles. Bhat- 
tal et al. used a constant e-softening with e ~500 AU, so that they could not resolve the 
formation of an object until its size became ^500 AU. 

Our simulation evolves with very small time-step due to the large number of particles. 
After about 10,000 time-steps, it becomes extremely inefficient to continue the simulation, as 
even if the time-step did not decrease any further, after another 10,000 time-steps the time 
would only progress by about 1% of the total time that far. Therefore, we have terminated 
our simulation at that point, after 0.64 Myr. There are ~260,000 particles at the end. 

The fact that our simulation has to stop very early prevents us from presenting figures 
with the time evolution of properties such as the mass, angular momentum and radius of 
fragments or the binding energies of groups of protostars. 

After ~0.61 Myr, a network of filaments forms within the shocked layer. Protostellar 
objects start forming at the intercepts of the filaments. On-the-fly splitting is applied at 
~0.63 Myr. In total, four objects have formed via fragmentation of the filaments when the 
simulation terminates. Two of them quickly merge. At the end of the simulation (~0.64 
Myr), three disc-like objects are still accreting matter from the filaments. They are already 
centrally condensed. They are rotating fast and are attended by thin spiral arms. An offset 
develops between the opposing accretion flows and the inflowing material steadily increases 
the speciflc angular momentum of the discs. They are spinning increasingly fast. At this 
early stage of the evolution there is no sign of interaction between the spiral arms and the 
accretion flows. We believe that secondary companions to the discs may eventually form in 
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our simulation via interaction between the spiral arms and the accretion flows. However, we 
note that the accretion flows in our simulation are extremely homogeneous, and not lumpy 
as they are in the simulation of Bhattal et al. 

The three discs are still in the sub-solar mass range, in accordance with the mass evolution 
of Fig. 3 of Bhattal et al. Specifically, the masses of the three disc-like objects are mi=0.35 
Mq, 7712=0.30 Mq and m3=0.19 Mq. They all extend out to ~120 AU. Their central densities 
have reached ff^"'^=2.1 x 10~^^ g cm~^, p^°'^=b.b x 10~^^ g cm~^ and p^^°'^=Q.S x 10~^^ g 
cm~^, which implies that the temperature at their centres is increasing. Their separation is 
of the order of 10^ AU and they appear to be a weakly bound system at this stage. The left 
panel of Fig. 15. II is a column density plot viewed along the z-axis and the right panel of Fig. 
15.11 is a column density plot viewed along the y-axis showing the three protostellar discs at 
the end of the simulation. 

Although it is possible that rotational instabilities may subsequently produce secondaries 
to the objects formed in our simulations, the initial formation of multiple protostars is due 
to layer fragmentation] this is not present in the simulation of Bhattal et al. The total mass 
of the three fragments of our simulation is of the same order of magnitude as the mass of the 
central object in the simulation of Bhattal et al., at corresponding times. 

Moreover, a network of 3-4 filaments is produced in the shocked layer of our simulation; 
Bhattal et al. do not produce a network. At the end of the simulation, the minimum density 
of the filaments is pfii=3.0 x 10"^'' g cm~^, which corresponds to n/^j ^5 x 10^ cm~^ (Fig. 

EH). 

We presume that the evolution of the density field is followed better in our simulation. 
The fact that the Jeans condition is obeyed, prevents a central object from forming artificially 
before the filaments. Since with on-the-fly splitting there is no preferred length scale, we can 
trust the detailed evolution of the discs themselves as well as their dynamical interaction. 
However, because of the large number of particles, the time-step decreases rapidly and we 
have to stop the simulations at a very early stage before the group of protostellar discs can 
evolve into a bound system. 

5.2.2 Simulation with Mo=75 Mq and b=OA 

For this simulation, Bhattal et al. found that a multiple system formed both via fragmentation 
of a single filament in the shocked layer and via rotational instability in the central, and more 
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The figure is provided separately 
due to its large size 



Figure 5.3. Column density plot for the 6=0.4 collision of two 75 M© clumps viewed along the z-axis 
at the end of the simulation, t -^0.69 Myr. The linear size of this plot is 0.008 pc. The colour table 
has units of 2.95 x 10^ g cm~^. 

massive, fragment formed in the filament. The fragments formed a bit later than the central 
object in their 6=0.2 collision. 

Their findings are repeated in our simulation with the exception of the rotational insta- 
bility event. 

As in the 6=0.2 collision, we terminate this simulation after ~0.69 Myr, since the time- 
step has decreased considerably and continuing becomes computationally inefficient. The 
number of particles has increased to 260,000 at the end of the simulation. 

A single filament forms perpendicular to the collision axis, at ~0.66 Myr. At the same 
time, on-the-fly splitting initiates. Four objects have formed in the filament by t ^^0.69 Myr 
and two more have started condensing out of the filament. The masses of the four objects 
are mi=0.33 M©, m2=m3=m4=0.02 M©. The more massive object is disc-like with a radius 
of ~170 AU. The three smaller ones are still spherical with radii of ~45 AU. Their central 
densities have reached /9?^"*^=2.0 x lO'^^ g ^^-3^ p|^"'==3.3 x 10-^^ g ^.^^-3^ p|^"'==l.3 x 
10-13 g cm-3 and ff^"'^=%.2 x 10-^^ g cm'^. 

The four objects are formed in random positions along the filament and are falling to- 
gether. They would presumably interact if the simulation were continued further. 



138 



CHAPTER 5. CLUMP-CLUMP COLLISIONS 



The figure is provided separately 
due to its large size 



Figure 5.4. Column density plot for the 6=0.5 collision of two 75 M0 clumps viewed along the z-axis 
at the end of the simulation, t ~0.75 Myr. The linear size of this plot is 0.02 pc. The colour table has 
units of 4.73 x 10'^ g cm"^. 

The more massive object rotates fast and might well become rotationally unstable at a 
later stage. It forms first and, after becoming self-gravitating, it evolves with the minimum 
time-step. This is why the other fragments seem to evolve rather slowly compared to it. 

Fig. 15.31 shows the four objects within the filament at the end of the simulation. At this 
stage, the minimum density of the filament is pjj/=2.0 x 10"^'' g cm~^, corresponding to 
71^2 ^5 X 10^ cm^'^. 

5.2.3 Simulation with Mo=75 M© and 6=0.5 

For this collision, Bhattal et al. found that a single filament formed diagonally in the x- 
y plane. The filament subsequently produced several fragments that merged. After ~1 
Myr, a binary with equal mass components was spiraling in towards a possible merger. The 
filament fragmented later than in their 6=0.4 collision. When they repeated their simulation 
with higher resolution, they resolved each binary component to be a binary itself, i.e. an 
hierarchical quadruple. 

With our simulation, we confirm their results. We stop the simulation at ~0.75 Myr, 
when it involves ~342,000 particles. The time-step decreases so much that it is inefficient to 
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continue any further. Particle splitting initiates at ~0.72 Myr. At ~0.73 Myr the filament 
forms. It soon starts fragmenting into several objects. Two of them, at the bottom right hand 
corner of Fig. 15.41 condense out first and eventually collapse faster, thereby dominating the 
smallest time-step bins. At ~0.746 Myr they merge. At the end of the simulation, another 
three objects appear to be forming (the two green lumps at the centre of Fig. 15.41 plus 
another one out of Fig. 15.41 above and to the left). 

All fragments remain spherical and rotate but considerably slower than in the two pre- 
vious simulations. In particular, they become very centrally condensed. Some other density 
enhancements appear in the filament suggesting that more fragments may be produced (via 
fragmentation of the filament). The masses of the four objects at the final stage are mi=0.46 
Mq, m2=0.07 Mq, m3=0.06 Mq and 772,4=0. 04 Mq. The most massive one has radius ~35 
AU and the radius of the other three is ~100 AU. Their central densities are p^^'*^=1.6 x 



^5 X 10^ cm ^. 
5.2.4 Conclusions 

Apart from the 6=0.2 simulation, our results are in general agreement with those of Bhattal 
et al. (|1998|) . In the 6=0.2 collision, we do not obtain a central object. Instead, a network 
of filaments forms, the filaments fragment, and a group of protostellar discs forms at the 
intercepts of the filaments. We believe that with our simulation, where fragmentation is 
properly modelled, we obtain a more realistic evolution for the shocked layer. 

We conclude that in all cases fragmentation of the filaments produces the seed for a group 
of protostars. Accretion-induced rotational instability can then produce secondaries to the 
protostars formed. Our simulations have to stop at an early stage due to the small time-step 
with which the simulation progresses; the time-step drops due to the large number of particles. 
Thus, we can not confirm the development of rotational instabilities. We expect that such 
instabilities may develop at a later stage when the angular momentum of the discs increases 
further. We anticipate that accretion-induced rotational instabilities are more probable with 
small 6 (6 <0.5). 

The filaments fragment later with increasing 6, as the leading edge of each colliding clump 
has to travel more into the other clump before the shocked layer becomes massive enough. 
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Bhattal et al found that the primary mass was larger with decreasing b. One might argue 
that this is evident in our simulations as well, since the total mass of the primaries is largest 
in the 6=0.2 collision (^^0.85 Mq). However, our evidence is somewhat inconclusive since 
only the fragments that collapse first are properly modelled; they dominate the small time 
bins and are evolved with the smallest time-steps while all other fragments are in suspended 
animation. 

The simulations of Bhattal et al. were of particularly low resolution and may have suffered 
from spurious fragmentation. Furthermore, we believe that since they could not properly 
model the density field (they could not resolve the filament formation in the 6=0.2 colli- 
sion), the primary masses they inferred were over-estimated. The primary protostars of our 
simulations are only a few ten thousand years old, since they form simultaneously with, or 
shortly after, the filaments. From our simulations a mass accretion rate > 10~^ Mq yr~^ can 
be inferred. With average mass accretion rates of 10^^ Mq yr~-^ for the first few hundred 
thousand years and 10"^ Mq yr~^ subsequently until they become a million years old, the 
primaries will end up as 3-6 Mq stars if they are not disturbed by rotational instability and/or 
mergers. If we consider as typical the 4-6 primary fragments that seem to be produced in 
our simulations, then 20-30 Mq in total will end up in stars in such collisions. This would 
imply that the efficiency of star formation in such clump-clump collisions is less than, or 
of the order of, 20% which is in accordance with the findings of Hunter et al. (1986) from 
simulations of colliding gas flows. In fact, they suggest a star formation efficiency slightly 
higher than this, but since their flows move faster than our clumps, their shocks are stronger, 
and thus fragmentation of the resulting filaments may produce more primaries. The star for- 
mation efficiency estimated from our simulations is similar to that of most molecular clouds 
( Rengarajan 198'5l ). The efficiency inferred from the simulations of Bhattal et al. is of the 
order 30-40% which is unrealistically high. 

Furthermore, we believe that the conclusion of Bhattal et al that in low-6 collisions 
unequal mass systems are produced whereas in medium- to high-6 collisions equal mass 
systems occur, is not necessarily sound since they arrive at this conclusion having not resolved 
the filament formation in low-6 collisions or the subsequent filament fragmentation. 

In all collisions, our simulations end before we can have conclusive evidence on the proto- 
stellar systems formed being bound or not. In all cases, the accretion flows bring them closer 
together steadily. 
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It is also interesting to note that in all cases the filaments produced have n/^j ^10^ cm^^, 
which means that such filaments could be observed in molecular line radiation using tracers 
that are excited at such densities, e.g. NH3, CS. 

Finally, we conclude that in low b collisions, a network of filaments develops instead of the 
single spindle-like filament formed in medium to high b collisions (6 ^0.4). In low b collisions, 
more mass from the two clouds ends up in the shocked layer, as they collide almost head-on. 
This increases the surface area of the layer. 

5.3 Low-mass clump collisions 

In our low-mass clump collision simulations, two equal mass clumps, with Afo=10 Mq, collide 
with equal but opposite velocities. We use Mq=W Mq, oq ~0.35 km s^^ (corresponding to 
ro=35 K); M=5, 10, 15 (corresponding to Ucw=0.5, 1.0 and 1.5 km s"^), 6=0.0, 0.2, 0.4, 
0.6 and 0.8, and 15,000 particles per clump. Table 15.21 gives a summary of all simulations 
conducted as well as the most important results. 

All collisions were conducted with and without application of particle splitting in order 
to be able to quantify the effect of particle splitting - i.e. the Jeans condition being obeyed 
- on the final result. All simulations evolve with very small time-step. Especially after 
the application of particle splitting, the time-step decreases due to the increased number of 
particles (decreased inter-particle distance; see Eqn. I2.46() . This happened in the higher mass 
clump collisions as well f ^5.2() . We stop the particle splitting simulations when it becomes 
inefficient for them to be continued any further (see discussion in ^5.2.1|) . We also stop the 
simulations without particle splitting at a time close to that of the corresponding particle 
splitting simulation only for comparison purposes. 

In the sequel all figures refer to the particle splitting simulations only. 

5.3.1 Simulation with Mo=10 Mq, M=5 and 6=0.0 

The initial conditions for the head-on collisions are shown in the top panel of Fig. 15.51 

This is a head-on collision with the clumps moving with relatively slow velocities {vcioud=0.5 

km s^^). On-the-fly splitting initiates at tspi ~0.455 Myr. The simulation stops at ~0.476 

Myr. At the end of the simulation there are 42,700 active particles. 

Only a single spherical fragment has formed. It forms at ~0.471 Myr. The fragment forms 
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Table 5.2. Summary of simulations and most important results for the low-clump collisions (Mo=10 Mq) for the different values of b and A4. 
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at the centre of the simulation and it accretes material from the shocked layer. No filament 
forms in this collision. At the end of the simulation, the fragment mass is 0.33 M© and its 
radius ~50 AU. Its central density has reached Ppeofe=4.75 x 10"^^ g cm~^. It rotates but not 
as fast as the fragments in the simulations of Bhattal et al. In this simulation, there seems to 
be (cylindrically) isotropic inflow of matter onto the fragment as opposed to accretion along 
filaments in the simulations of Bhattal et al. 

The simulation without on-the-fly splitting produces similar results. The protostar forms 
at the same time as above (~0.471 Myr). At the end of the simulation (~0.476 Myr), its 
mass is 0.47 M©. The Jeans condition is not obeyed above Pmax=^ x 10~^^ g cm~^. The 
protostar's radius is smaller (~32 AU) as it has become more centrally condensed (Ppeafc=l-75 
X 10~^^ g cm"^). We presume that with on-the-fly splitting we obtain a more realistic picture 
for this simulation even for such an early time in its evolution. 

5.3.2 Simulation with Mo=10 M©, M=10 and 6=0.0 

For this head-on collision the shock is stronger as the clumps move with higher velocities 
{vcioud=^-^ km s~^) than in the previous simulation. On-the-fly splitting initiates at tgpi 
~0.343 Myr. The simulation stops at ~0.370 Myr. At the end of the simulation there are 
45,500 active particles. 

Two filaments form at ^^0.360 Myr almost perpendicular to each other. At the intercept 
of the two filaments a single object forms and it accretes material from the two filaments. At 
the end of the simulation, the object is spherical and rotating. Its mass is 0.35 M© and its 
radius ~45 AU. Its central density has reached /9peafe=4.9 x 10"-*^^ g cm~^. 

The mass and the radius of the fragment are similar to those of the fragment of the 
previous simulation. The formation of the fragment happens faster in this simulation due to 
the fact that the two clumps are moving with higher velocities (this also forces the time-step 
of the simulation to decrease considerably at an earlier time). 

The simulation without on-the-fly splitting evolves in a similar way. In particular, two 
filaments form at ~0.360 Myr and the single object, formed on their intercept, is accreting 
material from them. It eventually becomes disc-like. At the end of the simulation (~0.367 
Myr), its mass is 0.27 M© and its radius ~38 AU. However, it has become more centrally 
condensed (/9peafc=7.1 x 10^^^ g cm~^). 
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5.3.3 Simulation with Mo=10 M©, M=15 and 6=0.0 

This is the head-on coUision involving the highest clump velocities (fdoMd=l-5 km s^^). On- 
the-fly splitting starts at tgpi ~0.310 Myr. The simulation stops at ~0.332 Myr. At the end 
of the simulation there are 44,100 active particles. 

A network of filaments forms at ~0.320 Myr and a single rotating object forms at the 
intercept of the filaments almost simultaneously with them. The central protostellar object is 
surrounded by an extended disc. Material from the filaments swirls onto the disc. Due to the 
high clump velocities and the fact that accretion to the object happens in a non-axisymmetric 
fashion, the angular momentum grows faster than in the two previous simulations. The 
angular momentum increases with time but, at the end of the simulation, it is not sufficient 
for the growth of rotational instabilities. 

At this time {t ~0.332 Myr - bottom panel of Fig. 15. 5|) . the mass of the protostellar 
object is 0.47 Mq and its radius ~145 AU. Its central density has reached Ppeak— 

2.9 X 10"^^ 

g cm~'^. The mass of the fragment is again similar to the masses of the protostars formed in 
the previous simulations. The fragment again forms faster than in the previous simulations, 
as the shocked layer forms faster due to the clumps colliding at a higher relative velocity. 

In the simulation without particle splitting, a network of filaments forms and the rotating 
protostellar disc-like object forms almost simultaneously with the filaments, at ~0.324 Myr. 
The protostellar disc rotates faster with time. After a short period of time, spiral arms 
appear in the disc. Up to the end of the simulation (~0.334 Myr), no significant interaction 
between the spiral arms and the accretion flows is observed. At this time, the mass of the 
protostar is 0.60 Mq, its radius ~100 AU, and it has become more centrally condensed 
(/jpea/c=l-2 X 10^^^ g cm^^). The protostar grows in mass considerably faster since the 
Jeans condition is not obeyed above Pmax=^ x 10^^^ g cm^^. Again, comparison of the 
the simulation without particle splitting and the on-the-fly splitting simulation demonstrates 
that the results obtained with particle splitting are different, and presumably more realistic. 

5.3.4 Effect of changing M. with constant 6=0.0 

Comparison between the three head-on collision simulations with different clump velocities, 
shows that protostars with similar mass form. The simulations stop at a very early stage of 
the accretion process, only a few thousand years after the protostars form. At these early 
stages of protostellar collapse, protostars appear to be accreting mass with an accretion rate 
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Figure 5.5. Column density plots for the 6=0.0 and A4=15 collision of two 10 M© clumps. Top: 
Initial conditions viewed along the z-axis. The linear size of this plot is 0.48 pc. The colour table has 
units of 8.20 x 10^ g cm~^. Bottom : Column density plot viewed along the x-axis at the end of the 
simulation, t ~0.332 Myr. The linear size of this plot is 0.035 pc. The colour table has units of 1.54 
X 10^ g cm~^. 
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Figure 5.6. Column density plots for the 6=0.2 and A4=5 collision of two 10 M© clumps. Top: 
Initial conditions viewed along the z-axis. The linear size of this plot is 0.56 pc. The colour table has 
units of 6.03 x 10'^ g cm~^. Bottom : Column density plot viewed along the z-axis at the end of the 
simulation, t ~0.496 Myr. The linear size of this plot is 0.016 pc. The colour table has units of 7.38 
X 10^ g cm~^. 
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of ~5 X 10~^ Mq yr^-*^. This value is similar to that found in the 75 Mq clump collisions 
f ^5.2|) . The protostars formed are accreting with a rate similar to that of Class objects. 

We also note that the angular momentum of the protostars increases with increasing 
clump velocity, as the protostar in the M=15 collision forms a well-defined disc, in contrast 
to the protostars in the other two simulations that remain spherical. There is no evidence 
for the growth of rotationally induced instabilities up to the time the collisions end. It is 
more likely for these kind of instabilities to occur in the simulation with Vcioud=^-^ km s~^, 
at a later time. Subsequent interaction of the spiral arms with the accretion flows may form 
companions to the primary protostar. 

The filaments become more well-defined with increasing clump velocity (i.e. as the shocks 
become stronger), and they start forming from a lower density. This implies that filaments 
are the mechanism by which material accretes onto the protostars in high compression shocks. 
We only present a column density plot for the simulation with Vcioud=^-^ km s^^ as this is the 
most interesting simulation having the densest filaments (with minimum density of pfii=2.5 
X 10~^^ g cm^^, i.e. ^5 x 10^ cm~^.). 

5.3.5 Simulation with Mo=10 M©, M=5 and 6=0.2 

The initial conditions for the 6=0.2 collisions are shown in the top panel of Fig. 15.61 

In the particle splitting simulation of the 6=0.2 and M=5 collision, on-the-fiy splitting 
starts at tspi ~0.463 Myr. A single filament forms perpendicular to the collision axis at 
~0.480 Myr. The filament feeds a single protostellar disc. The disc rotates fast. There is 
spiral structure in the disc. At the end of the simulation (~0.496 Myr), the spiral arms are 
very well-defined but there is no evidence for them interacting with the uniform accretion 
flows. We can not rule out the formation of secondaries at a later stage. 

The mass of the protostellar object is 0.70 Mq and its radius ~90 AU (bottom panel 
of Fig. 15. 6j) . Its central density has reached Ppeak=T-^ x 10~^^ g cm~'^. At the end of the 
simulation there are 55,300 active particles. 

The simulation without particle splitting produces very similar results. The single rotating 
disc-like fragment is again within a single filament perpendicular to the collision axis (the 
filament forms at ~0.480 Myr). Spiral arms develop in the disc. The fragment is not evolved 
so much as in the particle splitting simulation since the simulation without particle splitting 
ends a few thousand years earlier than the particle splitting simulation (at ~0.490 Myr). The 
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mass of the fragment at the end of the simulation is 0.45 Mq and its radius is ~84 AU. Its 
central density has reached /9peafc=5-4 x 10^^^ g cm^'^. We note that, although there is no 
evidence for interaction between the accretion flows and the spiral arms, the accretion flows 
in the simulation without particle splitting are lumpy. This is the main difference with the 
particle splitting simulation where the accretion through the filament is very smooth and 
without any lumps. 

Comparison between the particle splitting simulations with 6=0.0 and 6=0.2 (both with 
M=5) shows that the single protostar rotates faster with larger 6; the angular momentum 
is larger due to the increased impact parameter of the collision. However, both simulations 
end rather soon after the protostar formation and this prevents the angular momentum from 
becoming very large. The simulation with 6=0.2 ends at a later time as it takes more time 
for the shocked layer to become massive enough to fragment. 

5.3.6 Simulation with Mo=10 Mq, M=10 and b=0.2 

In the particle splitting simulation of the 6=0.2 and A^=10 collision, on-the-fly splitting 
starts at t^pi ~0.339 Myr. A network of filaments forms at this time. At ~0.360 Myr, two 
objects form at the intercepts of the filaments. The protostars rotate fast and accretion discs 
form around them shortly after their formation. At the end of the simulation (~0.396 Myr), 
the two protostellar discs are attended by spiral arms. There is no evidence for interaction 
between the spiral arms and the accretion flows. The discs are approaching each other, 
towards a possible capture or merger. The separation at periastron and alignment of the two 
approaching protostars, are the major factors which will determine whether the two objects 
merge or are captured into a binary orbit. It is too early to determine the evolution of this 
binary. 

The masses of the two objects at the final stage are mi =0.59 Mq and m2=0.41 Mq (top 
panel of Fig. 15.7(1 . The more massive protostar has radius ~76 AU and the radius of the other 
object is ~103 AU. The former is more centrally condensed (i.e. its disc is smaller). Their 
central densities are p^°'^=4:A x 10~^^ g cm~^ and p^°'^ = l.S x lO"^'^ g cm~^. The minimum 
density of the filaments is pfii=2.8 x 10^^'' g cm~^, i.e. 71^2 ^5 x 10^ cm^'^ (bottom panel 
of Fig. 15. 7j) . There is a possible third object starting to form a few time-steps before the end 
of the simulation. There are 64,300 active particles at the end of the simulation. 

The simulation without particle splitting produces similar results but can be followed 
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Figure 5.7. Column density plots for the 6=0.2 and A4=10 collision of two 10 Mq clumps at the end 
of the simulation, t -^0.396 Myr. Top: Column density plot viewed along the z-axis. The linear size 
of this plot is 0.016 pc. The colour table has units of 7.38 x 10® g cm~^. Bottom : Column density 
plot viewed along the y-axis. The linear size of this plot is 0.02 pc. The colour table has units of 4.72 
X 10® g cm~^. 
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Figure 5.8. Column density plots for the 6=0.2 and A4=15 collision of two 10 Mq clumps viewed 
along the x-axis at the end of the simulation, t ~0.368 Myr. Top: The linear size of this plot is 0.032 

pc. The colour table has units of 1.85 x lO*" g cm~^. Bottom : Zooming on the protostar on the left 
of the top panel. The linear size of this plot is 0.008 pc. The colour table has units of 2.95 x lO"" g 
cm~^. 
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only up to ^^0.382 Myr. In particular, a network of filaments forms at ~0.339 Myr and at 
~0.360 Myr a single object forms within the filaments. It rotates fast and becomes disc-like 
followed by spiral arms. There is no evidence for the formation of secondaries in the disc. The 
protostar is evolving with tiny time-step. This is the reason for the final time being shorter 
than that of the particle splitting simulation. The mass of the protostar at the end of the 
simulation is 0.44 Mq and its radius is ~110 AU. Its central density has reached ppeafc=3.25 
X 10~^^ g cm^'^. Towards the end of the simulation, two more objects start condensing out. 
We presume that the results of the on-the-fly splitting simulation are more realistic on two 
counts: two objects are formed instead of one, and the time-step is not so small allowing the 
simulation to evolve a bit longer. 

Comparison between the particle splitting simulations with 6=0.0 and b = 0.2 (both with 
Ai=10) shows that both protostars in the 6=0.2 collision are more massive than the single 
protostar in the 6=0.0 collision, as they are evolved for a longer time. In the former simulation 
both protostars rotate faster than the single object of the latter simulation. It appears that 
the angular momenta of the protostars increase with increasing 6. The 6=0.2 simulation ends 
later than the 6=0.0 as it takes more time for the clumps to move through each other before 
the shocked layer becomes massive enough to fragment. 

In our 6=0.2 collision involving 75 Mq clumps, the Mach number is 9. Comparison of the 
6=0.2 collisions with different mass clumps, shows that in both cases filament fragmentation 
is the mechanism for the formation of the primaries. Fewer primaries and filaments form 
in the low- mass clump collision as the shocked layer has lower surface density, S^. is 
proportional to the time elapsed from the beginning of the collision and in the 75 Mq clump 
collision the filaments form long after the end of the 10 Mq clump collision. Higher surface 
density gives smaller Jeans length in the layer (Aj oc Ti~^, Whitworth et al. (|1994b|) ). The 
mean separation between filaments in the layer and fragments in the filaments is of the order 
of the Jeans length. This is why this mean separation is smaller in the 75 Mq clump collision 
and hence more filaments form in the layer and more fragments in the filaments. 

5.3.7 Simulation with Mo=10 Mq, M = lb and f>=0.2 

In the particle splitting simulation of the 6=0.2 and 7W=15 collision a network of filaments 
forms at ~0.306 Myr. The filaments are more well-defined than before and they form on 
the whole surface of the shocked layer, not just at its centre as in previous simulations. 



152 



CHAPTER 5. CLUMP-CLUMP COLLISIONS 



On-the-fly splitting starts at ~0.336 Myr. Two objects form within the central filaments at 
~0.348 Myr. The protostars rotate fast and accretion discs form around them shortly after 
their formation. At the end of the simulation (~0.368 Myr), the disc of the more massive 
protostar is attended by spiral arms. There are density enhancements at the points where the 
spiral arms interact with the accretion flows (bottom panel of Fig. 15. 8|) . The enhancements 
in density may subsequently form one or two companions to the primary. 

It is too early to determine whether the two primaries form a bound or an unbound 
system. The masses of the two objects at the final stage are mi=0.53 Mq and m2=0.42 
Mq (top panel of Fig. 15. 8|) . There is also another 0.1 Mq associated with the spiral arms 
that will be the upper limit for the mass of the secondaries forming via rotational instability. 
The more massive protostar has radius ~42 AU and the radius of the other object is ~35 
AU. Their central densities are p^°'^=1.2 x 10^^^ g cm~^ and f}^"''^=?>.l x 10"^^ g cm~^. 
The minimum density of the filaments is pfu ~10^^^ g cm^'^, i.e. ~10'^ cm~^. There 
is a possible third object starting to condense in another filament. There are 70,200 active 
particles at the end of the simulation. 

The simulation without particle splitting produces similar results. It is followed up to 
~0.360 Myr. In particular, a network of filaments forms at ~0.306 Myr and at ~0.348 Myr 
a single object forms within the filaments. It rotates fast and becomes disc-like with spiral 
arms. There is no evidence for interactions between the accretion flows and the spiral arms. 
However, the protostar rotates faster than all other objects formed in simulations without 
particle splitting. Again, the protostar is evolving with very small time-step. This is the 
reason for the final time being shorter than that of the particle splitting simulation. The 
mass of the fragment at the end of the simulation is 0.45 M© and its radius is ~92 AU. Its 
central density has reached /9peafe=2-2 x 10^^^ g cm^^. Towards the end of the simulation, 
two more objects start condensing out. Again, we presume that with on-the-fly splitting we 
obtain more realistic results. 

When we compare the simulations with 6=0.0 and 6=0.2 (both with 7W=15), we conclude 
that the simulation evolution is delayed, and the angular momentum of the discs is increased 
with larger 6. A similar conclusion was reached for A^=10. 
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Figure 5.9. Column density plots for the 6=0.4 and A4=5 collision of two 10 Mq clumps viewed 
along the z-axis at the end of the simulation, t -^0.557 Myr. Top: The linear size of this plot is 0.024 
pc. The colour table has units of 3.28 x 10® g cm~^. Bottom : Zooming on the protostar on the 
bottom right hand corner of the top panel. The linear size of this plot is 0.01 pc. The colour table 
has units of 1.89 x 10^ g cm~^. 
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Figure 5.10. Column density plots for the 6=0.4 and A4=W collision of two 10 M© clumps viewed 
along the z-axis at the end of the simulation, t ~0.507 Myr. Top: The linear size of this plot is 0.048 
pc. The colour table has units of 8.20 x 10® g cm~^. Bottom : Zooming on the protostar on the 
bottom right hand corner of the top panel. The linear size of this plot is 0.008 pc. The colour table 
has units of 2.95 x 10'' g cm~^. 
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5.3.8 Effect of changing M. with constant 6=0.2 

By comparing the three simulations with 6=0.2, we arrive at similar conclusions to those of 
^5.3.41 The protostars appear to be accreting mass with an accretion rate of ^^5 x 10~^ 
yr"^ 

We also note that the angular momenta of the protostars increases with increasing clump 
velocity, as the protostars form more well-defined discs and rotate faster. A secondary com- 
panion may be starting to form at the end of the Vcioud=^-5 simulation. 

The filaments increase in number and become more well-defined with increasing clump 
velocity. Since oc Vdoud and Aj oc the Jeans length in the layer decreases with 

increasing Mach number of the collision (see discussion in ^5.3.6|) . 

The accretion rate does not increase with increasing Vdoud^ which means that the filaments 
are somehow holding up temporarily material from falling on to the protostars to release it 
at a later stage. The densest filaments were formed in the 7W=15 simulation (with minimum 
density pju ~10~^^ g cm~^, i.e. uh^ ~10^ cm~^). 

5.3.9 Simulation with Mo=10 Mq, M=5 and 6=0.4 

The initial conditions for the 6=0.4 collisions are shown in the top panel of Fig. 15.111 

In the particle splitting simulation of the 6=0.4 and M=5 collision on-the-fiy splitting 
starts at ~0.525 Myr. A single filament forms diagonally on the x-y plane at ~0.532 Myr. 
Two objects form towards the ends of the filament at ~0.535 Myr. The protostars rotate 
fast and accretion discs form around them shortly after their formation. At the end of the 
simulation (~0.557 Myr), one of the two protostellar discs is attended by strong spiral arms 
(bottom panel of Fig. 15. 9() . There are density enhancements at the points where the spiral 
arms interact with the accretion fiows. A secondary companion to the primary protostar 
might subsequently form from these enhancements. 

The masses of the two objects at the final stage are mi=0.67 Mq and m2=0.44 Mq (top 
panel of Fig. 15. 9() . The more massive protostar has radius ~77 AU and the radius of the 
other object is ~52 AU. Their central densities are =5.25 X 10-12 

g cm -3 and ^^^^=7.94 

X IQ-i^ g cm~'^ (the second object is more centrally condensed). It can not be determined 
whether the protostellar system is bound or not; the separation is ~3500 AU. The minimum 
density of the filaments is pfii=2.A3 x 10^^'' g cm~'^, i.e. tih^ ^5 x 10^ cm^^. There are 
79,400 active particles at the end of the simulation. 
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The simulation without particle splitting produces different results. It is followed to 
~0.544 Myr. In particular, a single thin filament forms at ~0.528 Myr and at ~0.533 Myr 
the filament starts fragmenting. It produces 10 objects. Two of them merge. They are all 
disc-like and they rotate fast, but there is no evidence for rotational instabilities. Again, the 
protostars are evolving with very small time-step. The total mass of the fragments at the end 
of the simulation is 1.17 Mq. The objects have radii of ~30 AU. Their central densities have 
reached ppeak=^ -0-7.68 x 10~^^ g cm~^. We presume again that with on-the-fly splitting 
we obtain more realistic results, as fragmentation of the filament into so many objects is 
suspicious. 

5.3.10 Simulation with Mo=10 Mq, M=10 and 6=0.4 

In the particle splitting simulation of the 6=0.4 and A^=10 collision on-the-fly splitting starts 
at ~0.432 Myr. A single filament forms at ~0.452 Myr. The filament is not as thin as in 
the previous simulation. One object forms at ~0.485 Myr. The protostar rotates fast and 
an accretion disc forms around it shortly after its formation. At the end of the simulation 
(~0.507 Myr), the disc is attended by spiral arms. Due to the density enhancements at the 
points where the spiral arms interact with the accretion flows, we conclude that it is likely 
that secondaries to the protostar would be formed subsequently (bottom panel of Fig. I5.1U|) . 
The mass of the object at the final stage is 0.48 Mq and its radius ~48 AU. Its central density 
is /Opeafe=3.4 X 10~^^ g cm~^. The minimum density of the filaments is pfii=2 x 10~^^ g cm~'^, 
i.e. riH^ ^5 x 10^ cm^^. There is a possible second object starting to condense at the other 
end of the filament (top panel of Fig. I5.1fl|) . There are 59,300 active particles at the end of 
the simulation. 

The simulation without particle splitting produces different results, similar to those of 
the simulation without particle splitting for 6=0.4 and M=5. It is followed up to ~0.477 
Myr. A single filament forms at ~0.458 Myr and at ~0.466 Myr the filament fragments. It 
produces 5 objects. Two of them merge. They are all disc-like and they rotate fast. There 
is no evidence for rotational instabilities growing in these protostars. Again, the protostars 
are evolving with very small time-step. The total mass of the fragments at the end of the 
simulation is 0.54 Mq. The objects have radii of 25-35 AU. Their central densities have 
reached /9peafc=0. 45-2.0 x 10^^^ g cm^'^. We presume that with on-the-fly splitting we obtain 
more realistic results. 



5.3. LOW-MASS CLUMP COLLISIONS 



157 



The figure is provided separately 
due to its large size 



The figure is provided separately 
due to its large size 



Figure 5.11. Column density plots for the 6=0.4 and A4=15 collision of two 10 Mq clumps. Top: 
Initial conditions viewed along the z-axis. The linear size of this plot is 0.024 pc. The colour table 
has units of 3.28 x 10^ g cm~^. Bottom : Column density plot viewed along the z-axis at the end of 
the simulation, t ^0.453 Myr. The linear size of this plot is 0.032 pc. The colour table has units of 
1.85 X 10^ g cm-2. 
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Figure 5.12. Column density plots for the 6=0.6 and A4=5 collision of two 10 M© clumps. Top: 
Initial conditions viewed along the z-axis. The linear size of this plot is 0.74 pc. The colour table has 
units of 3.45 x 10'^ g cm~^. Bottom : Column density plot viewed along the z-axis at the end of the 
simulation, t ~0.7 Myr. The linear size of this plot is 0.12 pc. The colour table has units of 1.31 x 
10^ g cm~^. 
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5.3.11 Simulation with Mo=10 Mq, M = 15 and 6=0.4 

In the particle splitting simulation of the 5=0.4 and M=15 collision no filaments are formed. 
On-the-fly splitting starts at ~0.397 Myr. One object forms at the centre of the domain at 
~0.433 Myr. The protostar rotates fast and an accretion disc forms around it. At the end of 
the simulation (~0.453 Myr), the protostellar disc is attended by weak spiral arms. There is 
no evidence for interaction between the spiral arms and the accretion flows. The mass of the 
object at the final stage is 0.31 Mq (bottom panel of Fig. 15.111) . Its radius is ~132 AU and 
its central density /?peaA;=3.3 x 10~^^ g cm~^. There are 39,700 active particles at the end of 
the simulation. 

The simulation without particle splitting produces similar results. It was followed to 
~0.437 Myr. A single object forms at ~0.428 Myr at the centre of the domain. It rotates 
fast. The protostar is evolving with very small time-step. The mass of the fragment at the 
end of the simulation is 0.32 Mq and its radius is ~97 AU. Its central density has reached 
yOpeaA;=2.43 X 10~^^ g cm~'^. We presume that with on-the-fly splitting we obtain more realistic 
results. 

5.3.12 Effect of changing Ai with constant 6=0.4 

Increasing M in simulations with 6=0.4 does not produce stronger shocks, nor does it increase 
the angular momentum of the protostars. The colliding clumps fail to form a strong shock in 
the A^=15 collision. This conclusion contrasts with our previous conclusion that increasing 
Ai gives stronger shocks with higher global angular momentum when b ^0.2. It appears that 
the 6=0.4 collisions are the borderline, as with b higher than 0.4 we shall see that the clump 
collisions do not form strong shocks. 

5.3.13 Simulations with Mo=10 M© and 6=0.6 

The initial conditions for the 6=0.6 collisions are shown in the top panel of Fig. 15.121 

In the particle splitting simulation of the 6=0.6 and A4=5 collision on-the-fly splitting 
starts at ~0.644 Myr. The two clumps move a long way into each other before the density 
increases significantly. Two single well-separated objects form at ~0.678 Myr. No filaments 
are formed in this collision. The protostars rotate fast and accretion discs form around them 
shortly after their formation. At the end of the simulation (~0.701 Myr), one of the two 
protostellar discs has strong spiral arms (bottom right protostar in the bottom panel of Fig. 
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I5.12() . It is possible that secondaries to the primary protostar will form shortly. The system 
of the two protostars is unbound as they are at a large distance (> 30,000 AU) and moving 
apart rapidly. 

The masses of the two objects at the final stage are mi =0.40 Mq and m2=0.28 Mq 
(bottom panel of Fig. I5.12|) . The more massive protostar has radius ~53 AU and the radius 
of the other object is ~85 AU. Their central densities are p^°'^=ZA x 10^^^ g cm~'^ and 
^peafc_g g ^ 10^^^ g cm~'^. There are 50,000 active particles at the end of the simulation. 

The simulation without particle splitting produces only one protostar (apparently corre- 
sponding to the bottom right protostar in the simulation with particle splitting; see bottom 
panel of Fig. I5.12|) . It forms at ~'0.671 Myr. It rotates fast and becomes disc-like with spiral 
arms. There is no evidence for interaction between the spiral arms and the accretion flows. 
Again, the protostar is evolving with very small time-step. The mass of the fragment at the 
end of the simulation (~0.684 Myr) is 0.30 Mq and its radius is ~80 AU. Its central density 
has reached yOpeaA:=2.3 x 10~^^ g cm~^. 

In the 6=0.6 simulations with ^A=1Q and A^=15 the clumps do not interact significantly 
to induce fragmentation, and this is also the case in all the 6=0.8 simulations f ^5. 3. 14(1 . 

5.3.14 Simulations with Mo=10 M© and 6=0.8 

The initial conditions for the 6=0.8 collisions are shown in the top panel of Fig. 15.131 There 
is no shocked layer formed in any of these collisions. The overlap of the two clumps is very 
small. We have conducted these simulations only for completeness. Here we present the 
end of the simulation with A4=5 and 6=0.8 (bottom panel of Fig. 15.131 t ~1.4 Myr). We 
have chosen to present this simulation because it is the one with the slowest moving clouds 
and, hence, the highest a priori probability that the collision could create a shock, due to 
gravitational focussing (as in H5. 3.13(1 . Particle splitting was not applied as the density in the 
two clumps did not exceed Pmax=Q x 10"^'^ g cm~^. 

5.4 Discussion 

We have conducted a series of cloud-cloud collision simulations using particle splitting. Two 
different sets of initial conditions were used ( ^5.1|) . In the set used previously by Bhattal et 
al, the mass of the colliding clumps was 75 Mq f ^5.2|) . In the second set of initial conditions 
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Figure 5.13. Column density plots for the 6=0.8 and A4=5 collision of two 10 Mq clumps. Top: 
Initial conditions viewed along the z-axis. The linear size of this plot is 0.72 pc. The colour table has 
units of 3.63 x 10^ g cm~^. Bottom : Column density plot viewed along the z-axis at the end of the 
simulation, t ~1.4 Myr. The linear size of this plot is 0.56 pc. The colour table has units of 6.03 x 
10^ g cm~^. 
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Clump Interaction | 


Jvl ]y Angular momentum J 
5 JJ- Filament Numbers | 
10 JJ- Filament Density | 

15 -IJ- tfrag i 


Transition 


No shock 



Table 5.3. Dependence of different quantities and phenomena on the increasing values of h and 
for simulations of cloud-cloud collisions. Note that the parameter space is divided in two sections: low 
b collisions produce stronger shocks. Large b collisions reduce the cloud interaction. The transition 
happens at 6=0.4. 

the clump mass was 10 f H5.3|) . 

With both sets of initial conditions, we have found that simulations with b ^0.5 produce 
shocked layers. Filaments or spindles form in the layers, with densities n/^j ^ 10^ cm~'^. 
Groups of protostellar discs form by fragmentation of the filaments. We have identified frag- 
mentation of the filaments as the common mechanism for Star Formation in these collisions. 
We have also found that the filaments act as the accretion channel that feeds the protostars 
with material. 

All protostars formed, have mass accretion rates of ~ 5 x 10~^ yr^^ for the first 10-20 
thousand years of their evolution. Their ages and mass accretion rates are comparable to 
those of Class objects. 

The simulations involving 75 Mq clumps, indicate that a group of 4-6 stars with masses 
between 3 and 6 forms as a result of these collisions. The inferred Star Formation 
efficiency is ~15-20%. To derive the star masses and the Star Formation efficiency, we assume 
the observed mass accretion rates of protostars and a total pre-main-sequence lifetime of about 
a million years. 

The 10 Mq clumps have smaller radii than the 75 M© clumps. Therefore, in 10 M© clump 
collisions, all the mass of the clumps enters the shocked layer in a shorter time than in 75 
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Mq clump collisions. Thus, accretion to the protostars through the filaments finishes earlier. 
Specifically, in low-mass clump collisions, 1-2 stars would form. Our simulations have 
shown that 1-2 protostars form in each collision. Thus, the Star Formation efficiency derived 
from low-mass clump collisions is 10-15%. 

Bhattal et al have found that accretion-induced rotational instabilities in protostellar 
discs can produce secondary companions to the protostars. In some of our simulations spiral 
arms form in the discs. We found evidence suggesting that the spiral arms are interacting 
with the accretion flows. However, our simulations stop at an early stage of the disc evolution 
due to time-step constraints. Thus, we can not confirm that such interactions are efficient in 
forming companions to the protostars. Formation of secondaries by accretion-induced rota- 
tional instabilities and/or disc-disc interactions would increase the Star Formation efficiency. 

Simulations with b ^0.5 suggest that the protostars are falling together along the fila- 
ments. The protostars could form bound systems. On smaller scales, they could merge or 
form binary systems by capture. To investigate their dynamic evolution, we need to use 
sink particles ( ^4.4.11 1) to replace the protostars, in order to prevent small time-steps form 
occuring. We have repeated the simulation with Mq=10 M©, A4=W and 6=0.2 by replacing 
the most massive protostar formed by a sink. The evolution of this simulation was similar to 
that of ^5.3.61 as the other protostar soon became dense enough to require small time-steps. 
We should develop an automated algorithm to replace all protostellar objects with densities 
above a certain threshold (e.g. Bonnell et al. p997|l ). 

To investigate the evolution of the system of protostars, we would also need to regulate the 
shear viscosity, for instance by using the time-dependent formulation of Morris & Monaghan 
H1997|l and/or the Balsara (|1995|1 switch. 

The detailed evolution of each simulation indicated a dependence on the values of the 
clump mass, clump velocity (i.e. collision Mach number) and impact parameter. Table 15.31 
summarises the dependence on the impact parameter, 6, and the Mach number, A4. With =^>, 
we indicate the increase in the values of these parameters. | indicates an increasing quantity 
or a phenomenon that becomes stronger and more frequent. | marks a decreasing quantity 
or a phenomenon that becomes less frequent. 

In particular, we find that with larger b, the angular momentum of the discs and the time 
for filament fragmentation are also larger. This makes rotational instabilities more probable 
for collisions with larger b. 
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Increasing the Mach number of the cohision, produces more fragments and networks of 
filaments, as the Jeans length in the shocked layer decreases. The filaments become more well- 
defined as the shock becomes stronger, and the time for filament fragmentation decreases. 
The angular momentum of the discs increases also with increasing M, indicating that in 
high Mach number collisions, rotational instabilities are likely to be more effective for the 
formation of secondary companions. 

Increasing the clump mass increases the number of filaments in the layer and the number 
of fragments in the filaments. It also creates stronger shocks and more well-defined filaments. 
The reason for this is again the decrease in the Jeans length of the layer. 

With 6 > 0.5 the colliding clumps interact less. For 6 > 0.5 and M > 5 the collisions do 
not create shocks. There seems to be a geometrical constraint on the strength of the shocked 
layers produced by cloud collisions: for low b collisions, the shocks are strongest as more 
mass enters in them. Collisions with 5=0.4 appear to produce less strong shocks. Therefore, 
we have concluded that the transition occurs at 6=0.4. From Table 15. HI we conclude that 
collisions with 6=0.2 and A^=10 or 15 are the most efficient. 

With particle splitting the Jeans condition remains obeyed at all times. The higher resolu- 
tion achieved with particle splitting decreases the values of the time-step and the simulations 
had to end at an early time. 

We presume that the particle splitting simulations produce more realistic results than 
the simulations of Bhattal et al. for the 75 Mq clump collisions and our 10 Mq clump 
collisions without particle splitting. Our presumption is supported by the fact that when the 
simulation with 10 M©, 6=0.2 and Ai=W was repeated with fine particles from the beginning 
(i.e. with 195,000 particles per cloud) and without particle splitting, the results were in very 
good agreement at comparable times. Repeating the above particle splitting simulation with 
clouds of 15,000 particles but with particles initially taken from a lattice, has given again 
results very similar to those of the particle splitting simulation of ^5.3.61 The simulation 
with particles initially taken from a lattice starts with less particle noise and it produces 
more symmetric results, i.e fragments in symmetric positions. However, the fragments have 
similar masses, radii and separations with those formed in the particle splitting simulation 
of ^5.3.61 In the future, we aim to test carefully the convergence of our results. 

Finally, our simulations have shown the formation of dense filaments with > 10^ 
cm~^. We predict that filaments, and in some cases networks of filaments, could be observed 
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in sites of dynamical star formation. In such sites, more than 2 or 3 Class objects would 
form almost simultaneously within a radius of few thousand AU. The filaments could be 
observed in NH3 molecular line observations. With arc-second resolution, filaments of the 
sizes inferred by our simulations could be observed in SFR regions at distances < 1 kpc. 
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Conclusions 



In this thesis we have studied Star Formation triggered by low-mass cloud collisions by means 
of numerical simulations. Many previous simulations of cloud-cloud collisions did not obey 
the Jeans condition for fragmentation due to their lack of sufficient numerical resolution. To 
overcome this limitation without using very large numbers of particles, we have developed 
particle splitting. This is a method that increases the number of particles only in places 
where higher resolution is required in order to satisfy the Jeans condition. 

The method has been applied to cloud-cloud collision simulations using clouds of different 
masses and examining collisions with different impact parameters and different relative cloud 
velocities. 

In this chapter, we summarise our conclusions. We have divided them in two sections: 
those on the numerical method and those on the results of cloud-cloud collisions. 

6.1 SPH and Particle Splitting 

For the numerical simulations conducted in this thesis we have used Smoothed Particle Hy- 
drodynamics, and Tree-Code-Gravity for the calculation of the gravitational forces. We have 
used a the second order Runge-Kutta time-integration scheme and multiple particle time 
steps and the Monaghan prescription for artificial viscosity (a = /3 = 1). 

The numerical code has been extensively tested in previous work, producing good results 
to most tests. We have shown that a) it is able to model collapse and b) it prevents numerical 
perturbations from developing in simulations of equilibrium spheres. The code can also model 
adequately shocks created by 5 collisions of gas flows. 
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The formulation of artificial viscosity in our code produces a large effective shear viscosity. 
As a result we cannot trust the detailed evolution of the discs around the protostars formed in 
our simulations. Future improvements to the code should seek to regulate the shear viscosity, 
for instance by using the time-dependent formulation of Morris & Monaghan (|1997|) and/or 
the Balsara (|1995|) switch. 

In this thesis, we have developed an algorithm for particle splitting, which increases the 
number of particles in a simulation, but only in regions where the resolution is not sufficient 
to continue obeying the Jeans condition. With the new method, a coarse particle in such a 
region is replaced by 13 particles of smaller mass. Our algorithm puts the 13 fine particles 
on the vertices of an equilateral lattice (fee). This is a lattice with minimum interstitial 
volume. It is appropriate for the 3-dimensional problems we model in this thesis. The density 
above which the Jeans condition stops being resolved is inversely proportional to the particle 
mass squared. By decreasing the particle mass 13-fold, the density at which the simulation 
reaches its resolution limit increases by 13^. In all the fine simulations performed in this 
thesis, adiabatic heating switched on at densities lower than the critical density. Therefore, 
the Jeans condition was obeyed throughout these simulations. In principle, we could split 
particles repeatedly, and thereby increase the resolution indefinitely. 

Identifying the coarse particles that need to be split is achieved with two separate versions 
of the new method. In nested splitting we manually decide the volume and the position of 
the region where the Jeans condition stops being obeyed, and we split all particles in this 
region. Subsequently, all coarse particles entering this region are split on-the-fly. In on-the- 
fly splitting we calculate the density above which the Jeans condition will be violated and 
all particles that obtain densities higher than this are split on-the-fly. On-the-fly splitting is 
preferred as particles are not split unnecessarily nor do we have to stop the coarse simulation. 

Particle splitting was found to introduce errors into the calculation of particle smoothing 
lengths. We had to revise the method for calculating particle smoothing lengths, so that 
it complied with the mixing of different mass particles. In calculating particle smoothing 
lengths we now take into account the amount of mass contained in a smoothing kernel, and 
not the number of neighbours. 

Both versions of particle splitting have been tested and give good results. SPH with parti- 
cle splitting was found to model adequetly isothermal collapse. In simulations of equilibrium 
spheres the errors introduced by particle splitting were not sufficient for perturbations to be 
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produced and propagated . 

Application of particle splitting to simulations of rotating clouds with m=2 density pertur- 
bations gives results that agree with the results of high resolution finite difference simulations. 
SPH simulations without particle splitting required at least twice as many particles to repro- 
duce the same results. Therefore, the new method is both efficient and economic, in terms 
of computational cost. Specifically, for the particle splitting simulations wc have used only 
~ 40% of the memory and ~ 25 — 30% of the CPU used in the 600,000 simulation without 
particle splitting. The simulations of rotating clouds with m=2 density perturbations are 
computationally demanding as a large fraction of the total mass ends up in the protostars 
formed. We expect particle splitting to be even more efficient in problems where a large frac- 
tion of the total mass is not evolved in fragmentation. On-the-fiy splitting is more economic 
than nested splitting and it is the version we use in simulations of cloud-cloud collisions. 

However, simulations of rotating clouds with m=2 density perturbations which start with 
a small number of particles, have shown that particle splitting, in response to the imminent 
violation of the Jeans condition, is only a necessary, and not a sufficient, condition for the 
reliability and efficiency of a simulation. 

Particle splitting can be applied not just to cloud-cloud collision simulations but also 
to other SPH simulations that also suffer from insufficient resolution. In particular, particle 
splitting could be used in collapse simulations, disc fragmentation calculations, simulations of 
disc interactions, simulations of stellar winds, modelling of mass transfer discs in cataclysmic 
variables and accretion discs in super massive black holes. 

6.2 Cloud-cloud collisions 

We have applied on-the-fly splitting to low-mass clump collision simulations. With particle 
splitting artificial fragmentation has been eliminated. Two sets of initial conditions have been 
used in our simulations. The first set repeat previous simulations which were performed by 
Bhattal et al. without particle splitting; the parameters are Mo=75 Mq, M=9 and 6=0.2, 0.4 
and 0.5, with 110,000 particles per clump. The second set explore a new set of parameters; 
Mo=10 Mq; M=5, 10, 15; 6=0.0, 0.2, 0.4, 0.6 and 0.8, and 15,000 particles per clump. 

The simulations with b ^0.5 produce shocked layers. Filaments or spindles form in the 
shocked layers, with densities nn^ ^ 10^ cm~^. 
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Systems of protostars are produced by fragmentation of the filaments. We identify frag- 
mentation of filaments as the common mechanism for Star Formation in these collisions. 

Most of the protostars are surrounded by discs. Rotational instabilities in these discs 
may produce secondaries to the primary protostars. Spiral arms are formed in almost all 
discs, but interaction between the spiral arms and the accretion flows is observed only in a 
few cases. However, our simulations end early due to the decreased time-step. Spiral arms 
in discs become more frequent with increasing b and M. 

The protostars formed in the b ^0.5 simulations are falling together along the filaments. 
At the end of the simulations, they are still at large separations, and we can not estimate if 
the system they form is bound and/or if they will merge at a later time. 

All protostars formed show mass accretion rates of 5 x 10~^ Mq yr~^ for the first 10-20 
thousand years of their evolution. Their ages and mass accretion rates arc comparable to 
those of Class objects. The inferred Star Formation efficiency is ^^15-20% for the 75 Mq 
clump collisions and ^^10-15% for the 10 Mq clump collisions. The difference arises from 
the number of protostars formed which is larger in the former collisions. Disc instabilities 
and disc-disc interactions can create low-mass companions to the primaries formed in our 
simulations, and thus slightly increase the inferred values of Star Formation efficiency. 

In high Mach number collisions {M. ^ 9) with b < 0.4, a network of filaments forms. The 
filaments have higher column densities in the 75 M© clump collisions. 

We predict that filaments, and in some cases networks of filaments, could be observed in 
sites of dynamical star formation. In such sites, more than 2 or 3 Class objects would form 
almost simultaneously within a radius of few thousand AU. The filaments could be observed 
in NH3 molecular line radiation. With arc-second resolution, filaments of the sizes inferred 
by our simulations could be observed in SFR regions at distances < 1 kpc. 

Due to time-step constraints, our simulations can only be followed for a few thousand 
years after protostar formation. Sink particles could replace the protostars as soon as they 
reach a certain density. An automated algorithm should be implemented for this purpose. 
With such an algorithm, we would not need to follow the detailed evolution of the individual 
discs and therefore the simulations could be followed further. We would not be able to 
monitor the efficiency of accretion-induced rotational instabilities and disc-disc interactions 
as mechanisms for secondary formation, but we could obtain more specific values for the 
Star Formation efficiency of the filament fragmentation mechanism induced by cloud-cloud 
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collisions. We could also determine whether the protostellar systems formed are bound. 

To obtain the detailed evolution of a disc and its interaction with the accretion flows from 
the filaments, we need to continue the simulations for a much longer time. Parallelisation of 
the code and increasing speed of computation achieved by super computers may make this 
possible. 
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Appendix A 

Jeans criterion of stability 



The fact that signatures of cohapse f Hl.31 2') only appear at certain sites in the interstellar 
medium, indicates that, in general, interstellar gas is in a quasi-static state, where the self- 
gravity is balanced by hydrostatic pressure, turbulence and possibly magnetic fields. We 
would like to know the point where this balance breaks, as at this point, collapse initiates 
and the gas is no longer efficiently supported. 

Jeans (|19()2|> dealt with the simple case of an infinite homogeneous gas at rest, supported 
only by its pressure. He concentrated on the velocity of wave propagation of a small fluctu- 
ation in density, 5/3, 

8v 

p— = -V6P + pV5^, (A.l) 
at 



with 



V^S^ = -47rGp. (A.3) 
(taken from Chandrasekhar ^1939^ ^. For adiabatic gas {6P = c^dp), Eqn. lA. II becomes 

d 

p— V • V = -c^V^5p + pV^5^, (A.4) 
ot 

or 

—dp = c^vHp + ATTGp6p. (A.5) 
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APPENDIX A. JEANS CRITERION OF STABILITY 



This is a typical wave equation, with solution of the form 

5p oc e'('^-^+"*\ 

as long as o"^ = c^k^ — AnGp. The velocity of propagation, or Jeans velocity, Vj, is therefore 



Vj = T = c[l 



AttGp\ 



1/2 



k \ c^k^ J 
The gas is unstable for all wave numbers k < kj, where 



kj = -UnGp)^/^. (A.6) 
c 

If we assume that the corresponding minimum wavelength 

corresponds to the radius of the smallest unstable spherical fragment, the mass of this frag- 
ment is 



^3^5/2 



which is called the Jeans mass^. 



similar formulation is given in Jeans II1929II . It is based on the excess of gravitational energy when 
collapse initiates. The arithmetic coefficients in both cases are of the same order of magnitude. 
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